Prepared for submission to JHEP 



JIMWLK evolution in the Gaussian approximation 



E. lancu'^ and D.N. Triantafyllopoulos^ 

"■Institut de Physique Theorique de Saclay, F-91191 Gif-sur-Yvette, France 
^ECT*, European Center for Theoretical Studies in Nuclear Physics and Related Areas, 
Strada delle Taharelle 286, 1-38123 Villazzano (TN), Italy 

E-mail: edmond.iancuOcea.fr, trianta@ectstar.eu 

Abstract: We demonstrate that the Bahtsky- JIMWLK equations describing the high-energy 
evolution of the n-point functions of the Wilson lines (the QCD scattering amplitudes in the 
r~| ! eikonal approximation) admit a controlled mean field approximation of the Gaussian type, for 

any value of the number of colors Nc- This approximation is strictly correct in the weak scatter- 
ing regime at relatively large transverse momenta, where it reproduces the BFKL dynamics, and 
in the strong scattering regime deeply at saturation, where it properly describes the evolution of 
the scattering amplitudes towards the respective black disk limits. The approximation scheme is 
fully specified by giving the 2-point function (the S-matrix for a color dipole), which in turn can 
be related to the solution to the Balitsky-Kovchegov equation, including at finite Nc- Any higher 
n-point function with n > 4 can be computed in terms of the dipole S'-matrix by solving a closed 
system of evolution equations (a simplified version of the respective Balitsky- JIMWLK equa- 
tions) which are local in the transverse coordinates. For simple configurations of the projectile 
in the transverse plane, our new results for the 4-point and the 6-point functions coincide with 
the high-energy extrapolations of the respective results in the McLerran-Venugopalan model. 
One cornerstone of our construction is a symmetry property of the JIMWLK evolution, that 
we notice here for the first time: the fact that, with increasing energy, a hadron is expanding 
its longitudinal support symmetrically around the light-cone. This corresponds to invariance 
under time reversal for the scattering amplitudes. 
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1 Introduction 

The final state of an ultrarelativistic hadron-hadron collision, as currently explored at RHIC and 
the LHC, is characterized by an extreme complexity in terms of the number and the distribution 
of the produced particles. The study of multiparticle correlations represents an essential tool 
for organizing this complexity and extracting physical information out of it. In particular, a 
recent measurement at RHIC of di-hadron correlations in deuteron-gold collisions [1] revealed 
an interesting phenomenon — the azimuthal correlations are rapidly suppressed when increasing 
the rapidity towards the fragmentation region of the deuteron — , which is qualitatively [2-5] 
and even semi-quantitatively [6, 7] consistent with the physical picture of gluon saturation in 
the nuclear wavefunction. For this interpretation to be firmly established, one needs a more 
precise understanding of the multi-particle correlations in the high-energy scattering and, in 
particular, of their evolution with increasing rapidity. This triggered new theoretical studies [8- 
12] of many-body correlations in the color glass condensate (CGC), which is the QCD effective 
theory for high-energy evolution and gluon saturation, to leading logarithmic accuracy at least. 

The central ingredient in the CGC theory is the JIMWLK (Jalilian-Marian, lancu, McLer- 
ran, Weigert, Leonidov, Kovner) equation [13-20], a functional renormalization group equation 
of the Fokker-Planck type which describes the non-linear evolution of the gluon distribution in 
the hadron wavefunction with increasing rapidity, or decreasing the gluon longitudinal momen- 
tum fraction x. When applied to an asymmetric, 'dilute-dense', scattering (like a pA collision). 
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the JIMWLK evolution can be equivalently reformulated as an infinite hierarchy of ordinary 
evolution equations, originally derived by Balitsky [21], which refer to gauge-invariant corre- 
lations built with products of Wilson lines. A 'Wilson line' is a path-ordered exponential of 
the color field in the target. It describes the scattering between a parton from the projectile 
(the proton) and the dense gluonic system in the target (the nucleus), in the eikonal approxi- 
mation. Via the optical theorem, the n-point functions of the Wilson lines can be related to 
cross-sections for particle production in pA collisions. For instance, the single-inclusive quark 
(or gluon) production is related to the S'-matrix of a 'color dipole' (the 2-point function of the 
Wilson lines). Similarly, the production of a pair of partons with similar rapidities is related 
to the 'color quadrupole' (the 4-point function). The suppression of azimuthal di-hadron cor- 
relations in d-|-Au collisions at RHIC [1] occurs in the right range of transverse momenta, of 
the order of the nuclear saturation momentum ~ 1 GeV, to be interpreted as a result of 
gluon saturation and multiple interactions in the scattering of the quadrupole. Such non-linear 
phenomena are mean field effects, which are likely to be correctly described by the JIMWLK 
evolution, although the latter is known to miss another important class of correlations — those 
associated with gluon number fluctuations in the dilute regime, or 'Pomeron loops'^ [23-27]. 

Motivated by the above considerations, there were several recent studies of the quadrupole 
evolution in the framework of the Balitsky- JIMWLK equations [10-12]. The results in Ref. [11] 
appeared as particularly intriguing. In that paper, one has numerically solved the JIMWLK 
equation by using its representation as a functional Langevin process [28] and used the results 
to evaluate the quadrupole 5-matrix for different rapidities and for special configurations of the 
4 external points in the transverse plane. Remarkably, the results thus obtained show a very 
good agreement with the heuristic extrapolation to high energy of the corresponding results in 
the McLerran-Venugopalan (MV) model [29, 30]. We recall that the MV model refers to a large 
nucleus {A ^ 1) at not too high energy (where the effects of the evolution are still negligible) 
and that in this model the CGC weight function is taken to be a Gaussian: the only non- 
trivial correlation of the color fields in the nucleus is their 2-point function, the 'unintegrated 
gluon distribution'. The 'high-energy extrapolation' alluded to above refers to using the MV 
expression for the quadrupole S-matrix in terms of the dipole S'-matrix [3, 8], but with the 
latter taken from the numerical solution to the JIMWLK equation at the rapidity of interest. 

Such extrapolations have often been used for phenomenological studies [3, 8, 31-35], but 
their justification from the viewpoint of the high-energy evolution remained obscure. A Gaussian 
Ansatz has also been used for mean field studies of the Balitsky- JIMWLK evolution [36-39] . But 
these previous studies have not convincingly addressed the issue of the validity of the Gaussian 
approximation — in particular, they did not justify its suitability for describing higher n-point 
functions, such as the quadrupole. (The only, qualitative, attempt in that sense is the 'random 
phase approximation' proposed in Ref. [39]; see the discussion in Sect. 3.4 below.) In principle, 
there is no contradiction between having a Gaussian weight function for the target field and 
still generating non-trivial correlations in the scattering of many-body projectiles: indeed, the 
scattering amplitudes are built with Wilson lines, which are non-linear in the target field to 

^The 'Pomeron loops' are formally higher order effects and, moreover, they are supressed by the running of the 
coupling — at least in the calculation of the single inclusive particle production [22]. However, there is currently 
no reliable estimate of their effects on correlations in multi-particle production. 
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all orders. But within the context of the JIMWLK evolution, such Gaussian approximations 
seem to be prohibited by the highly non-linear structure of the evolution equation, which is the 
mathematical expression of gluon saturation. 

In spite of this theoretical prejudice, the numerical results in Ref. [11] suggest that a Gaus- 
sian approximation to the JIMWLK evolution may nevertheless work. Another piece of evidence 
in that sense emerges from the recent analytic study in Ref. [12]. There, we have constructed 
an approximate version of the Balitsky-JIMWLK hierarchy which is simple enough to allow for 
explicit solutions. Then we have showed that, for the special configurations of the quadrupole 
considered in Ref. [11], these approximate solutions coincide with the respective predictions of 
the MV model extrapolated to high energy. But in that context too, the similarity with the 
MV model appears as merely an 'accident', with no deep motivation: the simplified hierarchy 
proposed in Ref. [12] is generated by the 'virtual' piece of the JIMWLK Hamiltonian, which is 
non-linear in the target field and therefore seems incompatible with a Gaussian approximation. 
Moreover, the approximations in Ref. [12] have been justified only in the limit where the number 
of colors Nc is large (formally, Nc — ?• oo). This does not explain the observation in Ref. [11] 
that the numerical solutions to the JIMWLK evolution for Nc = 3 are better reproduced by the 
finite-A'^c version of the MV model (with Nc = 3, of course) than by its large-A^^ limit. 

Our purpose in the present analysis is to clarify such 'coincidences' and 'apparent contra- 
dictions' by resolving the aforementioned tensions between the simplified hierarchy proposed 
in Ref. [12], the Gaussian approximation, and the large- Ac limit. The results that we shall 
obtain can be summarized as follows. We shall demonstrate that the JIMWLK equation admits 
indeed an approximate Gaussian solution for the CGC weight function, that this solution is 
unique within the limits of its accuracy, and that it is tantamount to a simplified system of 
evolution equations, which are linear (while being consistent with unitarity) and local in the 
transverse coordinates. In the limit where Ac — >• oo, these new equations reduce to those previ- 
ously proposed in Ref. [12]. The ultimate outcome of our analysis is a global approximation to 
the Balitsky-JIMWLK hierarchy, which is valid for any Ac and allows one to construct explicit, 
analytic, solutions for all the re-point functions of the Wilson lines. These approximate solu- 
tions are strictly correct in the limiting regimes at very large {k± ^ Qs{Y)) and, respectively, 
very small {kj_ <^ Qs(Y)) transverse momenta, and provide a smooth (infinitely differentiable) 
interpolation between these limits. Here, Qs{Y) denotes the saturation momentum in the target 
at a rapidity Y equal with the rapidity separation between the target and the projectile. 

To describe our results in more detail, let us first explain the distinction between 'real' 
and 'virtual' terms in the Balitsky-JIMWLK equations. The 'real' terms describe the evolution 
of the projectile via the emission of small-x gluons, whereas the 'virtual' terms express the 
probability for the projectile not to evolve, i.e. not to radiate such (small-x) gluons. The 
'virtual' terms dominate the evolution in the approach towards the unitarity (or 'black disk') 
limit, since in that regime the scattering is strong and the projectile has more chances to survive 
unscattered if it remains 'simple' — i.e., if it does not evolve by emitting more gluons. By 
the same token, the 'virtual' terms control the evolution of the many-body correlations which, 
within the context of JIMWLK, are built exclusively via non-linear effects (multiple scattering 
and gluon recombination) in the regime of strong scattering. More precisely, the 'real' terms are 
important for that process too — they include the non-linear effects responsible for unitarity 
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and saturation — , but deeply at saturation their role becomes very simple: they merely prohibit 
the emission of new gluons with low transverse momenta k± < Qs{Y)- Thus, one can follow 
the evolution of correlations at saturation by keeping only the 'virtual' terms in the Balitsky- 
JIMWLK equations, but supplementing them with a phase-space cutoff which expresses the 
effect of the 'real' terms. (This is strictly correct in a 'leading-logarithmic approximation' to be 
detailed in Sect. 3.4.) Moreover, since the simplified equations thus obtained are linear, they can 
be extended to also cover the BFKL evolution in the weak scattering regime at k_\_ S> QsiY)- 
Indeed, in that regime and to the accuracy of interest, the re-point functions of the Wilson lines 
reduce to linear combinations of the dipole scattering amplitude, with the latter obeying the 
BFKL equation. The BFKL dynamics involves both 'real' and 'virtual' terms, but it can be 
effectively taken into account by tuning the kernel in the 'virtual' terms — namely, by requiring 
this kernel to approach the solution to the BFKL equation at large k±. 

The above considerations, to be substantiated by the subsequent analysis, explain why it is 
possible to approximate the Balitsky-JIMWLK equations by simpler equations which are linear 
and whose overall structure is inherited from the 'virtual' terms in the original equations. Similar 
considerations have underlined our previous construction in Ref. [12], but their generalization 
to finite N^. (that we shall provide in this paper) turns out to be highly non-trivial. 

Another subtle aspect of our present analysis is the recognition of the fact that the simpli- 
fied equations that we shall propose (for either finite or infinite Nc) correspond to a Gaussian 
approximation for the CGC weight function. A priori, the association of a linear system of 
equations with a Gaussian approximation may look natural, but in the present case this is com- 
plicated by the fact that, as alluded to before, the 'virtual' piece of the JIMWLK Hamiltonian 
is non-linear in the target field to all orders. Such a non-linear structure seems to preclude 
any Gaussian solution. The resolution of this mathematical puzzle turns out to be interesting 
on physical grounds, as it sheds new light on the physical picture of the JIMWLK evolution. 
Namely, we shall show that the Wilson lines within the 'virtual' Hamiltonian do not represent 
genuine non-linear effects associated with saturation, rather they express the physical fact that, 
with increasing energy, the longitudinal support of the target expands symmetrically around the 
light-cone. That is, in contrast to a widespread opinion in the literature (see e.g. [28, 36- 
38, 40]), which was based on a misinterpretation of the mathematical structure of the JIMWLK 
equation, the gluon distribution in the target expands simultaneously towards larger and respec- 
tively smaller values of x~, in such a way to remain symmetric around x~ = 0. (We assume 
the target to propagate along the positive axis at nearly the speed of light and we define 
x~ = {x^ — x^)/^/2.) In turn, this symmetry has physical consequences for the multi-partonic 
scattering amplitudes: it implies that the re-point functions of the Wilson lines with re > 4 obey 
a special permutation symmetry — the mirror symmetry — which expresses their invariance 
under time reversal. 

To summarize, a Gaussian weight function which is symmetric in x~ and whose kernel is 
energy-dependent and interpolates between the solution to the BFKL equation at high trans- 
verse momenta k±_ ^ Qs and the JIMWLK (or 'dipole') kernel at low momenta k± <C Qs, 
provides a reasonable approximation to the JIMWLK equation, which is strictly correct in the 
limiting regimes alluded to above (for any value of Nc). Within its limits of validity, this approx- 
imation is essentially unique: different constructions for the kernel can differ from each other 
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only in the transition region around saturation, which is anyway not under control within the 
present approximation. 

In practice, it is convenient to trade this kernel for the dipole S'-matrix, which in turn 
can be obtained either as the solution to the Balitsky-Kovchegov (BK) equation [21, 41], or 
by solving a self-consistency condition similar to that in Ref. [36]. (The differences between 
these two expressions for the kernel should be viewed as an indicator of the stability of the 
approximation scheme.) Then the n-point functions with n > 4 (quadrupole, sextupole etc) can 
be determined in terms of the 2-point function (the dipole S-matrix) by solving the evolution 
equations associated with the Gaussian weight function. These equations become particularly 
simple at large Nc, where they reduce to the equations proposed in Ref. [12] and can be explicitly 
solved for arbitrary configurations of the n external points in the transverse plane. 

For finite Nc and for generic configurations, the equations are more complicated, as they 
couple the evolution of the various n-point functions with the same value of n. (For instance, 
the quadrupole mixes under the evolution with a system of two dipoles.) Yet, explicit solutions 
can be obtained under the simplifying assumption that the kernel of the Gaussian is a separable 
function of the rapidity and the transverse coordinates (plus an arbitrary function of Y; see 
Sect. 4.2 for details). This is certainly not the case for the actual kernel (say, as given by 
the solution to the BK equation), but it is a good piecewise approximation to it and it is 
furthermore true for the MV model^, that we shall take as our initial condition at low energy. 
So, not surprisingly, the expressions for the n~point functions that we shall obtain within this 
scenario are formally similar to the respective predictions of the MV model. One can reverse 
this last argument as follows: given that the Gaussian weight function is a good, piecewise 
approximation to the JIMWLK evolution, as we shall demonstrate, and that the kernel of this 
Gaussian can be taken to be separable within the relevant kinematical regimes, we expect the 
predictions of this approximation to be very close to those of the MV model extrapolated to 
high energy. 

For special configurations which are highly symmetric, exact solutions can be obtained at 
finite Nc even without assuming separability. We shall study various examples of this type for the 
4~point function and the 6-point function, and thus find some surprising factorization properties, 
that would be interesting to test against numerical solutions to the JIMWLK equation. For one 
particular configuration of the 6-point function, the exact numerical result is already known [11] 
and our respective analytic solution appears to agree with it quite well. 

2 The Balitsky— JIMWLK evolution equations 

In this section, we shall briefly review the general formulation of the JIMWLK evolution and 
then use the evolution equations satisfied by the dipole and the quadrupole 5-matrices in order 
to illustrate various properties of the evolution, which are important for what follows: the role 
and origin of the 'real' and 'virtual' terms, the factorization of multi-trace observables at large 
Nc, and, especially, the symmetric expansion of the longitudinal support of the target and the 
ensuing, 'mirror', symmetry of the n-point functions with n > 4. 

^By 'rapidity-dependence' within the MV model, we more precisely mean the dependence upon the longitudinal 
coordinate x~ . Within the JIMWLK evolution, there is a one-to-one correspondence between Y and x~ (see the 
discussion in Sect. 2.3 below). 
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2.1 JIMWLK evolution: a brief reminder 

The color glass condensate is an effective theory for the small-x part of the wavefunction of 
an energetic hadron: the gluons carrying a small fraction x ^ 1 of the hadron's longitudinal 
momentum are described as a random distribution of classical color fields generated by sources 
with larger momentum fractions x' ^ x. Given the high-energy kinematics, in particular the fact 
that the distribution of the color sources is 'frozen' by Lorentz time dilation, this color field can be 
chosen (in a suitable gauge) to have a single non-zero component, namely J^{x) = 6^~^aa{x~ ,x) 
for a hadron moving along the positive z axis^ (a 'right mover'). All the correlations of this 
field are encoded into a functional probability distribution, the 'CGC weight function' Vl^Y[a], 
which contains information about the evolution of the color sources with increasing 'rapidity' 

Y = ln(l/x), from some initial value Yq up to the value Y of interest. In the high energy regime 
where as{Y — Yq) > 1 and to leading logarithmic accuracy with respect to the large logarithm 

Y — Yq = ln{xQ/x), this evolution is described by a functional renormalization group equation 
for M^y[a], known as the JIMWLK equation. The latter can be given a Hamiltonian form, 

-^Wria] = HWY[a] , (2.1) 

where H is the JIMWLK Hamiltonian — a second-order, functional differential operator, whose 
most convenient form for the present purposes is that given in [42] and reads 

where we use the notation = J d^u. . . to simplify writing, Ai is the 'dipole kernel', 



{u — 1))^ 
{u — zy{z — v) 



and and V are Wilson lines in the adjoint representation: 



= Pexp 



ig J dx aa{x ,x)T° 



(2.4) 



with P denoting path-ordering in x~ . The above form of the Hamiltonian is valid only when act- 
ing on gauge-invariant functionals of Oa, which will always be the case throughout our analysis. 
In fact, the observables of interest are gauge-invariant products of Wilson lines (see below). 

The functional derivatives in Eq. (2.2) are understood to act at the largest value of x~ , that 
is, at the upper end point of path-ordered exponentials like that in Eq. (2.4) (see e.g. Eq. (2.11)). 
These derivatives do not commute with each other, but their commutator is proportional to 5uv 
(cf. Eq. (2.27)) and thus vanishes when multiplied by Aiuvz] hence, there is no ambiguity 
concerning the ordering of the functional derivatives in Eq. (2.2). One can also notice that 
the last two terms in the JIMWLK Hamiltonian, i.e. those proportional to VuVz and to VzVv 
respectively, are in fact identical to each other, as it can be checked by exchanging u ■v^ v and 
a <r^ b and by using the property Vz"''^ = for color matrices in the adjoint representation. 



^We use standard definitions for the light-cone coordinates: x'^ = {x'^,x ,x), with x^ = {t ± z)l\pi and 
X = {x,y). The field Ua is independent of the light-cone time a"*", because of Lorentz time dilation. 
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To fully specify the problem, one also needs an initial condition for Eq. (2.1) at y = Yq; at least 
for a sufficiently large nucleus, this initial condition is provided by the McLerran-Venugopalan 
(MV) model [29, 30] (see Sect. 3.2 below). 

Physical observables, like scattering amplitudes for external projectiles, are represented by 
gauge invariant operators 0[a] built with the field Oa, whose target expectation values are 
computed via functional averaging with the CGC weight function: 

{d)Y = j VaO[a]WY[a]. (2.5) 

By taking a derivative in this equation with respect to Y , using Eq. (2.1), and integrating by 
parts within the functional integral over a, one obtains an evolution equation for the observable, 
in which the JIMWLK Hamiltonian acts on the operator 0[a\ : 

'-^ = my. (2.6) 

Unlike the JIMWLK equation (2.1), this is not a functional equation anymore, but an integro- 
differential equation. However, due to the non-linear structure of the Hamiltonian (2.2) with 
respect to the field Eq. (2.6) is generally not a closed equation — the action of on O gen- 
erates additional operators in the right hand side — , but just a member of an infinite hierarchy 
of coupled equations — the Balitsky- JIMWLK equations. Although mathematically equiva- 
lent, the functional equation (2.1) and the Balitsky- JIMWLK hierarchy offer complementary 
perspectives over the high-energy evolution. Eq. (2.1) depicts the evolution of the target via the 
emission of an additional gluon with rapidity between Y and Y + dy, in the background of the 
color field a generated via previous emissions, at rapidities Y' < Y. The Wilson lines within 
the structure of the Hamiltonian (2.2) describe the scattering between this new gluon and the 
background field, in the eikonal approximation. The Balitsky- JIMWLK hierarchy rather refers 
to the evolution of the projectile, more precisely, of the operator which describes its scattering 
off the target. This scattering is again computed in the eikonal approximation, so the operator 
O is naturally built with Wilson lines — one such a line for each parton within the projectile. 

2.2 Evolution equations for the dipole and the quadrupole 

To be more explicit, we shall consider two specific projectiles: a 'color dipole' made with a 
quark-antiquark {qq) pair and a 'color quadrupole' made with two qq pairs. In both cases, the 
overall color state of the partonic system is a color singlet. The 5-matrix operators describing 
the forward scattering of these projectiles off the CGC target read 

Sxix2 = ^xiX2 = ^^{^xi^X2) : (2-7) 

for the color dipole and, respectively, 

QxiX2X3X4 = SxjX2XsX4 = ^ tl'(KciKK2Kc3^4) ' (2-8) 

for the color quadrupole. In these equations, and V are Wilson lines similar to those in 
Eq. (2.4), but in the fundamental representation. The results that we shall obtain for these two 
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partonic systems will be easy to extend to projectiles made with n qq pairs, for which 



As we shall see, within the high-energy evolution, such single-trace operators mix with the 
multi-trace operators, of the form 

o= ^tr(y4Kc2--0^tr(yJ^^^2---)^tr(y;^^.2---)- (2.10) 

In order to construct evolution equations according to Eq. (2.6), we need the action of the 
functional derivatives w.r.t. on the Wilson lines. This reads (with 6xu = S^'^^x — u)) 

^ V^ = ig6xut^Vl -^Vx = -i95xuVxt\ (2.11) 



"'-'it ""u 

By using these rules within Eqs. (2.6) and (2.2), it is straightforward to derive the evolution 
equations satisfied by S'-matrices for the dipole and the quadrupole. The respective derivations 
can be found in the literature (see e.g. the Appendix Ref. [12]), but here we shall nevertheless 
indicate a few intermediate steps (on the example of the dipole evolution), to emphasize the 
origin of the various terms in the equations. To that aim, it is useful to view the JIMWLK 
Hamiltonian (2.2) as the sum of two pieces, H = ffvirt + -f^rcab where -ffvirt corresponds to 
the first two terms in Eq. (2.2) and the -ffroal corresponds to the last two terms there. This 
division between 'virtual' and 'real' terms refers to the evolution of the projectile (see the physical 
discussion after Eq. (2.16) below) and should not be confounded with the corresponding division 
for the evolution of the target [36, 40]. By acting with these Hamiltonian pieces on the dipole 
^-matrix, one finds (with a = asNc/n). 

H^ixtSxiX2 — I ^ i\jo ) / •^^x\X2zSx\X2- (2.12) 



27r V A^c. 

and respectively (recall that the last two terms in Eq. (2.2), which define i/rcab are actually 
identical with each other) 



■^^XlX2Z \ SxizSzX2 A79 ^XlX2 )) (2.13) 



2yj- / ■'^xiX2Zy~-xiz^zx2 



c 



where the second line follows after reexpressing the adjoint Wilson line in terms of fundamental 
ones, according to 

[V^y^t^ = V^H^ = Vh"- V. (2.14) 
and then using the Fierz identity 

tY{t''At°-B) = ^trAtrS- ^tr(A5). (2.15) 



- 8 - 



By adding together the above results, one sees that the terms proportional to that would 

be suppressed at large Nc, exactly cancel between 'real' and 'virtual' contributions, and we are 
left with 

= J ■MxiX2z{SxizSzX2 ~ ^XiX2)Y ■ (2.16) 

This equation has the following physical interpretation: the first term in the right hand side, 
which is quadratic in S and has been generated by the 'real' piece of the Hamiltonian, cf. 
Eq. (2.13), describes the splitting of the original dipole {xi, X2) into two new dipoles {xi, z) 
and (z, X2), which then scatter off the target. More precisely, the evolution step consists in the 
emission of a soft gluon, so the original dipole gets replaced by a quark-antiquark-gluon system 
which is manifest in the first line of Eq. (2.13), but in large- A'^c limit (to which refers the first 
term in the second line of Eq. (2.13)), this emission is equivalent to the dipole splitting alluded 
to above. As for the second term in Eq. (2.16), i.e. the negative term linear in S which has 
been produced by iJvirtj it describes the reduction in the probability that the dipole survive 
in its original state — that is, the probability for the dipole not to emit. In what follows, we 
shall often refer to the terms produced by i^virt (-f^rcal) as the 'virtual' ('real') terms, but one 
should keep in mind that not all such terms are actually visible in the evolution equation in 
their standard form in the literature (to be also used in this paper): some of these terms may 
have canceled between 'real' and 'virtual' contributions. 

A similar discussion applies to the evolution equation for the quadrupole, which reads 



d{QxiX2X3X4}Y 

dY ~ 4^ 



X1X2Z + M -M X2X4Z^{SxizQ ZX2X-SX4)y 

J z 

~l~ {,'J^X\X2Z ~l~ ■J^X2XSZ ■J^X\Xiz){S zX2QxiZXSX4}Y 

~\~ {,.f^X2XsZ ~l~ .^^X'^X4Z .f^X2X4z){Sx^zQx\X2ZX4lY 

~\~ {^J^X\X4Z ~\~ J^X-iXiZ J^X\Xj,z){SzX4Qx\X2X'!,z')Y 

i,-^^XlX2Z ~l~ J^X^X^Z ^ ■^^XlX4Z ~^ ■^^X2X3z){QxiX2X3X4)y 

(,.^^XlX2Z ~l~ .^^X'^X4Z '^XlX-^Z .^^X2X4z}{SxiX2^X-jX4)y 

(,-^^xiX4Z ~l~ ■^^X2Xsz -^^xix-^z .^^X2X4z'){Sx-iX2^xiX4)Y • (^•1'^) 

Namely, the terms involving {SQ)y in the right hand side are 'real' terms describing the splitting 
of the original quadrupole into a new quadrupole plus a dipole, and have been all generated by 
the action of the last two terms in the Hamiltonian (2.2). The 'virtual' terms involving {Q)y and 
{SS)y are necessary for probability conservation, and have been generated by the first two terms 
in the Hamiltonian. Once again, all the terms subleading at large Nc (as separately generated 
by -ffvirt and H^eai) have canceled in the final equation. 

The above features are generic: they apply to the evolution equations obeyed by all the 
single-trace observables like Eq. (2.9). As visible on Eqs. (2.16) and (2.17), these equations are 
generally not closed : they couple single-trace observables with the multi-trace ones. E.g., the 
equation for the quadrupole also involves the 4-point function {SS)y and the 6-point function 
{SQ)y, which in turn are coupled (via the respective evolution equations) to even higher-point 
correlators. The equations obeyed by the multi-trace observables exhibit an interesting new 
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feature: they involve genuine ^/N"^ corrections, as generated when the two functional derivatives 
in Eq. (2.2) act on Wilson lines which belong to different traces (see e.g. Appendix F in [43] for 
an example). At large N^, these corrections can be neglected and then it is easy to check that 
the hierarchy admits the factorized solution 

{d)Y^{^^tr{VlV^,...))^{^^tr{Vi^^^^^^ (2.18) 

provided this factorization is already satisfied by the initial conditions. Then the hierarchy 
drastically simplifies: it breaks into a set of equations which can be solved one after the other 
(at least in principle). Namely, Eq. (2.16) becomes a closed equation for {S)y (the BK equation 
[21, 41]), Eq. (2.17) becomes an inhomogeneous equation for {Q)y with coefficients which depend 
upon {S)y [3], and so on. In practice, however, the resolution of these equations is hindered by 
their strong non-locality in the transverse coordinates. So far, only the (numerical) solution to 
the BK equation has been explicitly constructed. 

2.3 The mirror symmetry 

In this subsection, we shall discuss a symmetry property of the JIMWLK equation, which has not 
been noticed in the previous literature and which has far-reaching consequences: the symmetry 
of the target field distribution (the CGC) under reflection in x~ . 

To start with, we shall identify a mirror symmetry in the evolution equation (2.17) for the 
quadrupole, that can be easily demonstrated in the large-A'^c limit, but is likely to hold for 
any Nc- (It does so, at least, in the Gaussian approximation that we shall later construct.) 
Specifically, if the quadrupole S-matrix {QxxX2xzx4)y is symmetric under the exchange of the 
two antiquark Wilson lines (that is, the Wilson lines at X2 and x^) at the initial rapidity Yq — a 
condition which is indeed satisfied within the MV model [8] — , then this symmetry is preserved 
by the evolution. That is, for any Y > Yq, one has 

{QxiX2X^Xi)Y — {QxiX4Xj,X2 )y- (2.19) 

A similar property holds for the exchange of the quark Wilson lines at xi and x-^, but this is not 
independent of Eq. (2.19), since ^^313x2^^10^4 — QxiX4,xzx2 the cyclic symmetry of the trace. 
To demonstrate Eq. (2.19), let us consider the respective an^z-symmetric piece: 

Qxl^2x,x, = ^ Hvl^V^^ykV^,) - ^V^^Vx^VlV^,)] . (2.20) 

By using Eq. (2.17), it is easy to see that the associated expectation value obeys the following 
evolution equation: 



di(cixiX2X'iX4)Y 

dY ~ 4^ 



j {■MxiX2Z ~^ -MxiXiZ ■^X2X4z){SxizQzx2X:jX4)Y 

+ {■MxiX2Z + ■Mx2X3Z ~ -^X^Xsz) {S ZX2Q XiZX3X4)Y 
+ {J^X2XiZ + J^XiXiZ ~ ■^X2X4,z){Sx3zQx^X2ZX4)Y 

+ {-J^xiXiZ + -J^x-iXiZ ~ ■^x\xs,z) {S zxiQ xiX2Xsz)y 

- iMx,X2Z + Mx,X4Z + Mx,X4Z + Mx2Xsz){Qx7^x,xJy] ■ (2.21) 
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At large Nc, where one can factorize (5(5''''^™)y ^ (5')y (Q''''^™)y , Eq. (2.21) becomes a ho- 
mogeneous equation which implies that {Q^^"^)y = at any Y provided this condition was 
satisfied at Yq. In turn, this implies the symmetry property (2.19). 

By inspection of the higher equations in the Balitsky-JIMWLK equations, one can check 
that a similar symmetry holds for all the n-point functions of the Wilson lines. For instance, 
the equation obeyed by the sextupole ^-matrix {S^^^)y is explicitly shown in Appendix B of 
Ref. [12]. From this equation, one can read the following symmetry property: 

)y • (2.22) 

The generalization of this property to the 2n-point function shown in Eq. (2.9) reads 

{^XlX2...X2n-2X2n-lX2„)Y = {SxiX2nX2„-lX2n-2---X2)Y ■ (2.23) 

To better understand the content of this symmetry, it is useful to give a pictorial representa- 
tion for it. To that aim, consider a generic configuration of the quadrupole operator QxiX2X3X4 in 
the transverse plane, as illustrated in Fig. 1, and join the four points by oriented lines, which fol- 
low the direction of color multiplication. In this way, one constructs a closed, oriented, contour, 
whose orientation indicates the flow of color within the operator. By repeating this procedure 
for the 'permuted' operator QxiX4X3X2: one obtains a similar contour, where however the orien- 
tation of the color flow is reversed. One can similarly check that, for a general n-point function, 
the symmetry property (2.23) refers to changing the contour orientation, say from clockwise to 
counterclockwise. Such a change would also result from the reflection in a mirror, so we shall 
refer to the symmetry property (2.23) as the 'mirror symmetry'. Additional arguments in the 
favor of this name will be given below. 

There are several reasons why this this symmetry is so important for us here. First, as 
we shall shortly argue, this corresponds to an important symmetry property of the scattering 
amplitudes: their invariance under time-reversal. Second, the way how this symmetry is actually 
preserved by the JIMWLK evolution is very interesting as it sheds light on the physical picture of 
the target field distribution: with increasing Y, the color glass condensate expands symmetrically 
around x~ = 0. Third, this symmetry will later guide us in the construction of a mean field 
approximation to the Balitsky-JIMWLK hierarchy. 

To understand the relation to time-reversal, let us present another pictorial representation 
for the two quadrupole operators which enter Eq. (2.19): this is shown in Fig. 2, where the 
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Figure 2. A different pictorial representation of the operators QxiXQ,xs,xi (left) and Qxix^x^x-i (I'ight), 
which emphasizes the fact that they get exchanged with each other via time reversal, with 'time' — . 



transverse space is schematically represented by the vertical axis, whereas the horizontal axis 
refers to x~ — the light-cone time for the projectile. The Wilson lines are now explicitly shown, 
as the oriented horizontal lines extending along the x~ axis and connected with each other, via 
matrix multiplication, at x~ — )• ±(X). Once again, the orientation of these lines corresponds to 
the direction of the color flow. Clearly, these two figures get exchanged with each other when 
inverting the arrow of time. In the left figure, corresponding to QxiX2X3Xi, the 4-body system 
starts at x~ — )■ —00 as a set of 2 dipoles, {xi,X2) and (2:3, 034), which then exchange color 
with each other at x" — )• 00 and thus reconnect into the new dipoles (0^1,0^4) and (333, a;2). 
In the right figure, the opposite process happens: the system starts with the dipoles [xi^x^) 
and (a;3,a;2)) which then reconnect at x~ — )• 00 into the dipoles (031,032) and (0:3, 034), thus 
yielding the quadrupole Qxixixsx2- Hence, the symmetry property (2.19) corresponds indeed to 
invariance under time-reversal, as anticipated. 

We now turn to the physical interpretation of the mirror symmetry in the context of 
the JIMWLK evolution. One can check that the symmetric structure of the virtual terms in 
Eq. (2.17) stems from the combined action of the first two terms in Eq. (2.2). Half of the 'virtual' 
terms are generated by the first term, proportional to the color unity matrix, but by themselves 
these terms do not show the mirror symmetry; this symmetry is recovered only after adding 
the other half of the 'virtual' terms, as generated by the second term in Eq. (2.2), proportional 
to Vilv^. As an example, consider two of the 'virtual' terms in the r.h.s. of Eq. (2.17) whose 
coefficients get exchanged with each other under the exchange 032 -H- 034 : Mxxxaz{Qxix2xzxa) 
and MxiX2z{QxiX2X3X4) ■ The first of them is generated when acting with the first term in the 
Hamiltonian on the pair (031,034) of the quadrupole, whereas the second one emerges from the 
action of the second term in H on the pair (031,032). 

Hence, to elucidate this symmetry, one needs to better understand the action of the JIMWLK 
Hamiltonian. As manifest from Eq. (2.11), the functional derivatives within the Hamiltonian 
act as generators of infinitesimal color rotations of the Wilson lines at their upper end point 
in x~] that is, they act as Lie derivatives for the color group SU(A''c). These color rotations 
express the evolution of the target color field aa(x~, 03) with increasing rapidity: performing one 
infinitesimal step in the evolution, from 1" to y + dy, amounts to 'integrating out' one layer of 
quantum fluctuations within the target wavefunction — the gluons with longitudinal momentum 
fractions between x = e~^ and x' = e~^^^'^^) — and results in adding one additional layer to 
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the classical color field aa{x~,x). The fact that the JIMWLK Hamiltonian acts on the Wilson 
lines via color rotations at the largest value in x~ means that the new layer of color field is 
located at larger x~ as compared to the previous layers. 

This argument makes it tempting to conclude that, with increasing Y, the support of the 
target field aa{x~ ,x) extends only towards increasing x~ , thus yielding a field distribution which 
is asymmetric in x~ . This was indeed the prevailing viewpoint in the original literature on the 
JIMWLK evolution (see e.g. [36, 40]), but now we shall argue that this is actually not quite 
right: although the functional derivatives in Eq. (2.2) have a one-sided action which amounts 
to color rotations at the largest value of x~ alone, the overall structure of the Hamiltonian is 
such that the target field is nevertheless built symmetrically in x~ . In fact, it is precisely this 
symmetry of the target field distribution under reflection in x~ which is responsible for the 
mirror symmetry in the evolution equations. 

To see that, it is useful to notice that Eq. (2.2) can be alternatively rewritten as [44, 45] 

2 

H = ^g^3 j -Muvzi^JtuJiv + JruJrv + "^^z^Jru-Jlv^ ^ (2.24) 

where and are functional differential operators acting as 'left' and 'right' Lie derivatives 
— that is, the generators of infinitesimal color rotations at the largest and, respectively, smallest 
value of x~. They are defined as 

\ S \ ^ S 

Jlu = ~~ T~;r 5 '^Ru = TT' ' (2.25) 

and satisfy 

Jlu vl = -d^u f'vl vl = s^^ v;t^^ (2.26) 

where the second equation follows from the first one after using Eq. (2.14). These equations 
imply the following commutation relations 

['^Lit' -^Lv] = ^f"'^'^ Jlu^uv , [Jru^ "^to] = ^/"'^'^Jru^uv, [Jlu^ ^Rv] = 0' (2.27) 

showing that the two sets of generators satisfy two independent SU(A''c) Lie algebras. 

The physical interpretation of the various terms in Eq. (2.24) is quite transparent: the 
action of on the quark Wilson line Vx (the first equation in Eq. (2.26)) amounts to the 
addition of an infinitesimal layer of color field at the largest values of x~ , whereas the action 
of J^^ is tantamount to a corresponding addition at the smallest values of x~ . Hence, the 
manifest symmetry of Eq. (2.24) under the exchange L ■(r^ R implies that, during the high- 
energy evolution, the distribution of the target color field a^{x~,x) — by which we mean its 
support and correlations — is built symmetrically in x~ around x~ = 0. Moreover, this is 
also the origin of the mirror symmetry since, as previously noticed, the latter follows from the 
combined action of the first two terms in the Hamiltonian (2.2), or (2.24) — those which get 
interchanged with each other under the permutation L o -R of the Lie derivatives. 

To better appreciate the differences between an evolution which is symmetric in x~ and 
one which is not, it is instructive to consider the evolution of a Wilson line, say for a quark 
projectile. When computing the target average (2.5) with the CGC weight function at rapidity 
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y, the support of the color field aa{x ,a;) is restricted to 
Xq exp(y — Yq). (This follows from the uncertainty principle: the softest gluon modes that have 
been integrated over have longitudinal momentum p'^ = xP with x = e~^, with P = the total 



hadron momentum; hence, they are delocalized in x over a distance 
quark Wilson line can be equivalently rewritten as 



Kf = P 



exp 



Iff 



dx aa{x ,x)t°' 



oc e^.) Thus, the 



(2.28) 



This makes it manifest that, with increasing y, the Wilson line 'grows' simultaneously at its 
both endpoints. To visualize the effect of one step in the evolution (y — )• y + dy), it is useful to 
discretize rapidity by writing = ne, with e an infinitesimal rapidity interval. Then under one 
additional step n — )• n + 1, the upper bound of the support extends as x~ = x^(y„) — )• x~j^i = 
x~(l + e) and the Wilson lines evolves as Vn — )■ V^+i, with 



Vn+ii^) = exp[ic/ean+i(a3)] y„^(a3) exp[ic/ea_(„+i) (a;)] 



(2.29) 



where 



an+i{x) = x^aa{x^j^^,x)t"- and a_(n+i){x) = x^aa{-x^^]^,x)t°' (2.30) 

represent the additional fields generated in this evolution step per unit of space-time rapidity. 
The infinitesimal gauge rotations associated with these new fields can be expanded in powers of 
e. Strictly speaking, this expansion must be pushed to quadratic order in e, to match with the 
fact that the evolution Hamiltonian (2.24) involves second order functional derivatives. However, 
the quadratic terms arising from the expansion of a given Wilson line do not contribute to the 
evolution of gauge-invariant observables: they would yield 'tadpole' contributions oc 6uv^ but 
the dipole kernel in the Hamiltonian vanishes when u = v. In other terms, the two functional 
derivatives within H must act on different Wilson lines within the observable to give a non-zero 
result. So, we can restrict the expansion of (2.29) to linear order in e, which yields 

Vl_,^{x) - V^{x) = ige[an+i{x) V,l{x) + V^{x)a_^^+,){x)] + 0{e^) . (2.31) 

Clearly, the two terms in the r.h.s. correspond to the infinitesimal, 'left' and 'right', color 
rotations in Eq. (2.26). If instead of the symmetric evolution above, one would have considered 
an asymmetric one, where the target fields expands towards positive values of x~ alone, the 
analog of Eqs. (2.29)-(2.31) would have involved the 'left' infinitesimal color precession alone. 

In Ref. [28], the JIMWLK evolution has been reformulated as a random walk in the space 
of Wilson lines, which is formally such that one additional step corresponds to an infinitesimal 
rotation of V^{x) on the 'left' alone. However, by inspection of the manipulations there, one 
can check that the additional contribution a„+i(a;) to the target field in the {n + l)th step is 
such that, in reality, that step simultaneously generate a color precession on the 'left' and on the 
'right'. That is, the Langevin process introduced in Ref. [28] does in fact describe a symmetric 
evolution for the Wilson lines (or for the target field distribution), although this has not been 
recognized there. 
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3 The Gaussian approximation 



In this section we shall demonstrate that the JIMWLK equation for the CGC weight function 
admits an approximate Gaussian solution which properly captures both the BFKL dynamics in 
the dilute regime at k± ^ Qs{Y) and the approach towards the black disk limit in the saturation 
regime at k± <C Qs{Y). Our analysis improves over previous, related, constructions in the 
literature [36-38] at two important levels: (i) we actually justify the Gaussian approximation 
— including for the description of the higher-point correlation functions and for finite Nc — , 
on the basis of the Balitsky-JIMWLK equations; (ii) we implement the 'mirror' symmetry 
discussed in Sect. 2.3, that is, we construct a Gaussian distribution which is symmetric in x~ at 
any Y. As we shall see, this last condition is in fact compulsory to achieve a faithful description 
of the JIMWLK dynamics deeply at saturation. 

The material of this section is organized as follows: the Gaussian weight function is intro- 
duced in Sect. 3.1, compared to the MV model in Sect. 3.2, and justified in Sects. 3.3 and 3.4 by 
comparison with piecewise approximations to the JIMWLK equations in the limiting regimes 
alluded to above. 

3.1 The Gaussian w^eight function 

The most general Gaussian weight function which is consistent with gauge symmetry"^ and 
describes a target field distribution which is symmetric in x~ reads 



Wy [a] = My exp 



-- / _ dx / aaix ,xi)'y {x ,xi,X2)aa{x ,3^2) 

^ J—x., JxiX2 



SY[a], (3.1) 



where x^(Y) = Xq exp(y— Iq) and the kernel j~^{x~ ,xi,X2) is an even function of x~ , assumed 
to be invertible. The functional (5-function 6y [a] ensures that the target field vanishes at larger 
longitudinal coordinates > xJj{Y) : 

SY[a]^ n nn'^(«-(^"'^))- m 

Here, 6{aa{x~ , x)) is the usual (^-function and a discretization of the space-time is understood. 
Finally, the overall normalization factor My in Eq. (3.1) is such that / VaWY[a\ = 1. 

Eq. (3.1) implies that the only non-trivial correlation of the target fields is their 2-point 
function, which is moreover local in x~ : 

{aa{x^ ,xi)ab{x2 ,X2))y = S"'^ Q{xJ.j{Y) - \x^\) 6{x:[ - X2)^ix^ ,xi,X2) . (3.3) 

Within this Gaussian approximation, the locality in x~ is required by gauge symmetry: to 
preserve the latter, any non-locality in x" should be accompanied by gauge links (Wilson lines) 
built with the field a", which would spoil Gaussianity. 

The Gaussian distribution (3.1) is manifestly symmetric in x~ around x~ = and this 
symmetry is preserved by the high energy evolution. In fact, Eq. (3.1) depends upon Y only via 



"^By 'gauge symmetry' we more precisely have in mind here the class of gauges within which the target color 
field has the structure A'^ = S'^'^Oa- Some gauge artifacts, which are inherent in Eq. (3.1) but turn out to be 
harmless in practice, will be later discussed. 
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the two endpoints, x^.j{Y) and —xJ,j{Y), of the support in x~ , meaning that the high-energy 
evolution proceeds via the symmetric expansion of the color field distribution towards both 
larger and smaller values of x~ . Specifically, by using the methods in Refs. [28, 36], one can 
check that the Gaussian weight function (3.1) obeys the following evolution equation: 



dWrja] 
dY 



I I 'yY{u,v) 

^ Juv 



+ 



Wvia] 



(3.4) 



where the 'left' ('right') functional derivatives act on the target field at the largest (smallest) 



(Y) and respectively x 



(Y). Also ^y{xi,X2) denotes 



value of X , that is, at x 
the field correlator per unit space-time rapidity as produced in the last step of the evolution, 



7y(a;i,a;2) = x,. 'y{±x..,xi,X2) . 



(3.5) 



Eq. (3.4) makes it manifest that the momentum rapidity and the space-time rapidity are iden- 
tified with each other by the high-energy evolution. In order to solve Eq. (3.4), one also needs 
the generalization of Eq. (3.5) to intermediate values y <Y for the space-time rapidity, that is 



7(X ,Xi,X2), 



with 



y = Yo+ln- 



(3.6) 



Eq. (3.4) should be viewed as a mean field approximation to the JIMWLK equation (2.1). 
It shows the same 'left -right' symmetry as the original equation, cf. Eq. (2.24), and hence it is 
consistent with the mirror symmetry discussed in Sect. 2.3. Clearly, this would not be the case 
if, instead of the symmetric Gaussian (3.1), one would consider an asymmetric one, say with 
support at < < xJ^{Y), as in the previous literature [36-38]: the corresponding evolution 
equation would contain only the 'left' functional derivatives — i.e., only the first term inside the 
brackets in Eq. (3.4). 

To justify the Gaussian Ansatz (3.1) for the CGC weight function, we shall shortly compare 
the associated evolution equation (3.4) to the actual JIMWLK equation, in different kinematical 
regimes. In this process, we shall deduce piecewise approximations for the kernel 7y(a;i,a;2), 
valid at high {k± » Qs{Y)) and low {k± <^ Qs{Y)) momenta, respectively. 

3.2 The McLerran Venugopalan model 

Before we turn to the JIMWLK evolution, let us briefly discuss the McLerran- Venugopalan 
(MV) model [29, 30] that we shall take as our initial condition at rapidity Yq. Besides providing 
the initial conditions, this model (and its ad-hoc extrapolation towards high-energy) will serve 
as a baseline of comparison for the mean-field results that we shall later obtain. Its discussion 
will also give us an opportunity to clarify some subtle aspects of the Gaussian approximation, 
like the gauge artifacts in the 'a-representation', cf. Eq. (3.1). 

In the MV model one assumes that the color charges in the nucleus are uncorrelated valence 
quarks. Accordingly the distribution of the color charge density p^{x~,x) is a Gaussian with a 
kernel which is local in the transverse plane: 



Wyo[p] = J^Yo exp 



dx 



p"'{x ,x)p"'{x ,x) 
X{x~,x) 



^Yq [p] , 



(3.7) 
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where of course X(—x ,x) = X{x ,x). Eq. (3.7) implies: 

{paix^ ,Xi)pb{x^ ,X2))y = 5abQ{x^ -\X^\) S{X^ - X2)Sx-,x^\{x^,Xi) . (3.8) 

The quantity X{x~ ,x) has the meaning of color charge squared per unit transverse area per unit 
longitudinal distance. In general the nucleus is assumed to be homogeneous in the transverse 
plane, i.e. the kernel in Eq. (3.7) is taken to be independent of x. Under that assumption, the 
calculation of expectation values in the MV model is not sensitive to the detailed dependence 
of the kernel upon x~ , but only to its integral 

_dx-\{x-) = 2 / dy\y= (3.9) 

which physically represents the color charge squared per unit area. In the above equation, the 
quantity \y = \x~\\{x~) (the strength of the charge correlator per unit space-time rapidity) has 
been defined by analogy with Eq. (3.6). The last equality follows after counting the color charges 
of the valence quarks within a nucleus with atomic number A and transverse area t^R\ (see e.g. 
[40]). The fact that it is only the integrated quantity (3.9) which matters arises from the fact 
that, under the present assumptions, the charge correlator (3.8) is separable as a function of x~ 
and the transverse coordinates. We shall return to this issue in Sects. 4.2 and 4.3. 

Eq. (3.7) is gauge invariant, but in order to make contact with the a-representation that we 
use throughout this paper, we shall henceforward consider it within the class of gauges where 
the target field is of the form Aa = d^^aa- Then aa{x^ x) is related to the color charge density 
Pa{x~ ,x) via the 2-dimensional Coulomb equation: —V\aa = Pa- So, for a homogeneous 
target, Eq. (3.8) implies the following expression for the 2-point function for the color field, in 
transverse momentum space (we denote r = xi — X2 and k± = \k\): 

%{k) = j d2re'^-7y(r-) = ^ for y < Yo . (3.10) 

Here and from now on, we prefer to work with the expressions of the various correlators per unit 
space-time rapidity., cf. Eq. (3.6), since these are the expressions which most directly enter the 
mean-field evolution equations like (3.4). 

Eq. (3.10) raises a potential problem: the Fourier transform of this expression back to the 
transverse coordinate space is not well defined, as it involves a (quadratic) infrared divergence. 
This problem reflects the fact that, by itself, the field aa(x~, x) is not invariant under the resid- 
ual gauge transformations which preserve the structure Aa = 5^^~^aa for the target field. The 
infinitesimal version of such a transformation reads aa{x~,x) — ?• aa{x~ ,x) + d'^uja{x~)., with 
Ua{x~) an arbitrary function [42]; so, clearly, the color charge density pa{x~ ,x) is invariant un- 
der this transformation. Strictly speaking, the general weight function Wy (and, in particular, 
its Gaussian approximation, Eq. (3.1)) should be written as a functional of /3, to make gauge 
symmetry manifest. On the other hand, observables like scattering amplitudes are built with 
Wilson lines, which are path-ordered exponentials of a. Taken separately, one Wilson line is 
not gauge invariant (rather, it transforms via color rotations [42]), but the physically relevant 
operators, which involve a product of such lines, cf. Eq. (2.9), are invariant. Whenever com- 
puting the expectation value of such a gauge-invariant operator, there is no problem with using 
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the weight function in the a-representation, as given in Eq. (3.1): all the gauge artifacts cancel 
out in the final result. 

As an example, consider the calculation of the dipole S'-matrix within the MV model. The 
corresponding result is well known and reads (see also Sect. 4.2) 

where we have assumed the MV model to apply at all the rapidities Y <Yo and we defined 

^Yo{xi,X2) = g^Cp I dy [%{xi,xi) + %{X2,X2) -2%{xi,X2)], (3.12) 

J —oo 

with Cp = (iV^ — l)/2Nc- Eq. (3.12) involves only the following linear combination of the target 
field correlators 

^y{x,,X2) ^ -%{X,,X2) + 1[%{X,,X,) +%{X2,X2)] = J ^[^ - ' ^"'^] ' (3-13) 

which is gauge-invariant, since under a residual gauge transformation the target field a"(x~, x) 
changes by a a;-independent quantity. The sign in the r.h.s. of Eq. (3.13) is such that ^y{xi, X2) 
be positive-semidefinite. The last equality in Eq. (3.13), which involves the color charge corre- 
lator Ay, illustrates the fact that the infrared divergences due to gauge artifacts cancel out in 
the linear combination (3.13). Strictly speaking, the above integral over k still has a logarithmic 
infrared divergence, but this is milder than the quadratic divergence appearing in the Fourier 
transform of %{k) in Eq. (3.10). The remaining divergence is not a gauge artifact anymore, 
but a 'physical' singularity of this model: it reflects the lack of correlations among the color 
sources. After taking into account the high-energy evolution, transverse correlations get built 
which screen out this divergence, as we shall shortly see. For completeness, let us estimate the 
final integral in Eq. (3.12): introducing an infrared cutoff A to regularize the remaining infrared 
divergence and writing r = \xi — X2\, one finds 

ry,{r)=g'Cp:^ / ,,,,2 .4 - ^In— ^, (3.14) 



2TrRl J (27r)2 kj_ 4 r^A'^ 

where Qq = 2a'^CFA/ R'j^ is essentially the nuclear saturation scale^ (as probed by a quark- 
antiquark dipole) at the initial rapidity Yq. Although obtained within the MV model, the above 
results are generic in the following sense: all the gauge-invariant observables computed in the 
Gaussian approximation involve the kernel jy{xi,xi) of the Gaussian (the correlator of the 
target color field) only via the linear combination shown in Eq. (3.13). So, in practice, there is 
no problem with using the a-representation, as shown in Eq. (3.1). 

Let us conclude this subsection with a remark on the calculation of expectation values 
within the MV model. The similarity between the respective weight function, Eq. (3.7), and 
the Gaussian approximation to the JIMWLK evolution, Eq. (3.1), makes it clear that one 
can consider the MV model as the result of a fictitious 'evolution' in which the target charge 
distribution is built in layers of x~ , from = up to = Xq . Specifically, let Wx-[p] 
denote the generalization of Eq. (3.7) in which Xq is replaced by X~ and assume the nucleus 



^More precisely, QsiYo) is defined by the condition Tygir = 2/Qs{Yo)) ~ 0{1). 
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to be homogeneous in the transverse plane. Then Eq. (3.7) is the solution to the following, 
functional, evolution equation (compare to Eq. (3.4)) 



dx- 




5 



6 



•Lv 



+ 



5 



5 



Rv 



,a 



Wx-[P\ 



(3.15) 



integrated from X~ = up to X~ = . In this equation, and /O^ refer to the color charge 
densities at x~ = X~ and x~ = —X~, respectively. When applied to the evolution of the 
Wilson-line correlations, Eq. (3.15) amounts to constructing the Wilson lines via symmetric 
iterations, i.e. via infinitesimal color precessions which proceed simultaneously 'on the left' and 
'on the right', as shown in Eq. (2.29). However, within the context of the MV model, this 
symmetric iteration is merely a choice of a discretization prescription and any other choice is 
equally good. As a matter of fact, the common choice in the literature in this context (see e.g. 
[3, 8, 31, 32]) is to perform asymmetric iterations 'on the left' : 



where this time n refers to a discretization of the x axis. This procedure is tantamount to 
solving the following evolution equation 



from X = —Xq up to X = Xq . In practice, one often takes Xq — t- oo, since the results are 



The above discussion sheds more light on the role of the 'left-right' symmetry in the evo- 
lution equations. So long as the CGC weight function is given (like in the MV model) and the 
associated evolution equations are merely used as a convenient device to compute expectation 
values, the symmetric discretization in Eq. (2.29) is not compulsory and it might not even be 
the most convenient one in practice. However, for the JIMWLK equation and any (mean field) 
approximation to it, the symmetric iteration is the only one to be correct, since this is how 
the target field distribution gets actually built via quantum evolution: the 'outer' layers (those 
located at larger values of |a;~|) are constructed after the 'inner' ones (those at smaller |x~|), and 
the new correlations built in one step depend upon the color field produced in all the previous 
steps. Hence, it would make no sense to consider an asymmetric evolution, like Eq. (3.17), since 
this would violate causality in the domain of negative x~ . 

3.3 Weak scattering regime: the BFKL dynamics 

In what follows, we shall study the JIMWLK evolution in two limiting regimes — large transverse 
momenta k± ^ Qs{Y) in this section and relatively small momenta k± <^ Qs(Y) in the next 
subsection — with the purpose of showing that, in both regimes, the evolution is consistent 
with a mean field approximation of the type shown in Eq. (3.4). We recall that Qs{Y) is the 
saturation momentum in the target (in a frame in which the target carries most of the total 
rapidity separation Y) and it increases with Y very fast. For a multi-point correlation function 
like the quadrupole (2.8), the statement that the 'transverse momenta are much larger than 
Qs^ means that all the transverse separations rjj = \xi — Xj\ between the external points are 




(3.16) 




(3.17) 
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much smaller than l/QgiY). Similarly, by 'momenta much smaller than Qs\ we mean that 
Tij 3> 1/QsiY) for any pair {xi,Xj) of external points. Very asymmetric configurations, where 
some of the distances rij are much larger than l/Qs{Y) while the others are much smaller, are 
strictly speaking not covered by the present analysis and must be separately studied. We shall 
discuss some examples of that kind in Sect. 4.4 below. 

For high transverse momenta k± ^ Qs(Y), the gluon density in the target is low, meaning 
that the corresponding color field is weak: g j dx~a ^ 1. It is then possible to expand the 
Wilson lines to lowest non-trivial order in the field in their exponent, within both the JIMWLK 
Hamiltonian and the operators defining the observables. For an operator like the dipole S'-matrix 
Eq. (2.7), we need to push the expansion in ga up to the second order, since the linear terms 
vanish after averaging. Introducing the dipole T-matrix operator Tx^x-z = 1 — SxiX27 whose 
expectation value represents the corresponding scattering amplitude, this expansion yields 

(T.ia,2)y =^ -OV with al^ j dx-aa{x~,x). (3.18) 

The weak scattering regime corresponds to (T)y ^ 1. Note that Eq. (3.18) involves only the 
linear combination (3.13) of the target field correlators, in agreement with the discussion in 
Sect. 3.2. The similar expansion for the quadrupole S-matrix, Eq. (2.8), yields 

1 {,Qx-iX2X'i,X4)Y — {TxiX2 '^XlX^ ~l~ Tx^X4 ~^ Tx2X-J Tx2X4 ~l~ Tx.^^X4)y J (3.19) 

where it is understood that {T)y is evaluated according to Eq. (3.18). More generally, in this 
dilute regime, all the n~point functions of the type shown in Eqs. (2.9) or (2.10) reduce to 
linear combinations of dipole amplitudes. This already shows that a Gaussian approximation 
for the CGC weight function should be indeed possible, to the accuracy of interest. To identify 
this approximation, let us also consider the weak-field limit of the JIMWLK Hamiltonian. Its 
obtention is facilitated by observing that Eq. (2.2) can be rewritten as 

1 + V^V^ - - Wv = (l - Vi%) (l - Wv) ■ (3.20) 

The leading order terms in the dilute regime are then obtained by expanding the Wilson lines 
within each of the two parentheses above to linear order in ga. (This amounts to an expansion 
of the original structure in the l.h.s. of Eq. (3.20) up to quadratic order.) For instance, 

l-Vl% -ig{al-al)T\ (3.21) 

with as defined in Eq. (3.18). After also using {T"')hc = i/"'"^) one finds H ~ -ffBFKL with 

This Hamiltonian is supposed to act on operators which are themselves evaluated in the weak- 
scattering regime and hence are quadratic functions of the field a^, as illustrated in (3.18) and 
(3.19). Clearly, the only evolution equation of interest for us here is that obeyed by the dipole 
scattering amplitude (3.18). This is readily obtained as 

d{Tx^x2)Y _ a f - ^ , . 

QY — 2yi- / ^13'2Z X^'^iClZ ' -'^2332 '''X1X2/Y' yO.^O) 
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and is recognized as the BFKL equation [46-48], that is, the equation obtained after hnearizing 
Eq. (2.16) with respect to {T)y- By using its solution, one can compute any other n-point 
function of the Wilson lines, like Eq. (3.19), in this dilute regime. 

We now construct the Gaussian approximation which reproduces the BFKL equation. To 
that aim, we shall compare the mean-field equation for {Tx-^x,^)y generated by Eq. (3.4) with 
Eq. (3.23) and thus deduce an approximate expression for ^y{u,v) valid in this linear regime. 
Notice that the left and right functional derivatives yield identical results when acting on the 
field a% which is integrated over x~ . Hence, Eqs. (3.4) and (3.18) imply 



dY ,,,, AN, A J^^"' \ 6ab 5a% ^"^^ "'^^^ / Y 



MFA 

= 25"'' [7y(a;i, a^i) + 7y(iC2, X2) - 27y (aii, 2:2)] 
- 1 

= 5'^^7y(^i,^2), (3.24) 
with '-^y{xi,X2) defined as in Eq. (3.13). The last equation can be integrated to yield 

f-y 

{Txix2)y =2g^CF fY{xi,X2) with fYixi,X2) = / dyjy{xi,X2). (3.25) 

MFA 

It is easy to check that the same expression for (T)y would be obtained by directly evaluating 
the expectation value in Eq. (3.18) with the help of Eq. (3.3). But its above derivation via 
the mean-field equation of motion has the merit to emphasize that the evolution equations for 
gauge-invariant observables generated by the Gaussian approximation involve the well-behaved 
kernel 7y (aji, 2:2) in spite of the fact that the corresponding functional equation (3.4) features the 
(generally ill defined) kernel 7y(a;i,a;2). This property is generic: it holds beyond the present, 
BFKL, approximation. Thus, for all practical purposes one can replace 7y(M,i?) — ?• — 7y(w,i7) 
within Eq. (3.4). This replacement works in the same way as that of the original kernel in the 
JIMWLK equation [17, 19, 20] by the dipole kernel in Eq. (2.2): the new kernel is to be used 
only when acting on gauge— invariant observables and it has the property to vanish at u = v. 

Returning to the mean-field expression (3.25) for (T)y, this must be consistent with the 
BFKL equation (3.23). This is clearly the case provided the function fY{xi,X2) itself satisfies 
the BFKL equation: 

dfY{xi,X2) _ Q_ 

dY ~ 2tt 

The initial conditions for the above equations can be taken from the MV model, which yields 
(for r < 1/Qo) : (r(r))y, = 2g^CF fY,{r) = ^Yo{r), with Fy^ given in Eq. (3.14). 

The solution to Eq. (3.26) is by now well understood. Here we will just remind that the 
BFKL evolution introduces transverse correlations between the 'color sources' (radiated gluons) 
which ensure that the solution /y (2:1,332) becomes infrared finite after a rapidity evolution 
Y — Yq ^ 1/ a. In particular, in the window for 'extended geometric scaling' [49-51], which holds 
for transverse momenta relatively close to (but still larger than) the saturation momentum 
Qs{Y)-, one has^ Ay(fc) = fc^7y(fc) oc k^^^ '^"'^ with 7^ 0.63 (the 'BFKL anomalous dimension 



/ MxiX2z[fY{xi,z) + fY{z,X2) - fY{xi,X2)] . (3.26) 
J z 



^Notice that A;*7y(fe) = A:^7i (fe) since in momentum space the difference between 71- (fe) and 7Y-(fc) is 
proportional to 5^'^\k). 
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at saturation'). Then, clearly, the integral over k in Eq. (3.13) is well defined when computed 
within the BFKL approximation. 

To summarize, the mean-field equation (3.4), where it is understood that the kernel can be 
replaced as jy — ^ — 7y with the function 7y(it, i?) determined by Eq. (3.26), properly encodes 
the BFKL evolution of the dipole amplitude in the weak scattering regime. This conclusion 
holds for any value of the number of colors Nc and it extends to all the n-point functions like 
(2.9) and (2.10) which, in this regime, reduce to linear combinations of dipole amplitudes. 

Before concluding this section, let us recall that there are also other aspects of the BFKL 
dynamics, which cannot be encoded into a Gaussian weight function. They refer to operators 
more complicated than those in Eq. (2.9), which already at weak scattering involve more than two 
gluon exchanges; that is, to lowest order in the weak field expansion, they involve polynomials 
in a of a degree higher than two. (Such operators can be obtained e.g. by subtracting the 
dipolar contributions to the Wilson-line operators in Eqs. (2.9)~(2.10).) An example of that 
type is the 'odderon' operator, which describes C-odd exchanges and which in perturbation 
theory starts with three gluon exchanges. The corresponding evolution equation is correctly 
encoded (to leading logarithmic accuracy) in the JIMWLK equation [42] — in particular, its 
low-density limit, known as the 'BKP equation' [52-54], is generated by the weak-field limit 
(3.22) of the JIMWLK Hamiltonian [42] — but this description goes beyond the purpose of a 
Gaussian approximation, which by construction can encode only the 2-point function of the a 
field. For instance, to describe odderon effects in the initial conditions, one needs an extension 
of the MV model allowing for a non-trivial 3-point function [55]. 

What is however remarkable about the Gaussian approximation that we pursue here is its 
capacity to encode non-trivial correlations among n Wilson lines with arbitrary n in the strong 
scattering regime, where the linear relation between the n-point functions and the 2-point 
function does not hold anymore. This will be discussed in the next subsection. 

3.4 Strong scattering regime: the dominance of the 'virtual' terms 

For relatively low transverse momenta kj_ < Qs(Y), the gluon occupation numbers in the target 
wavefunction saturate at a large value of order l/a^, meaning that g J dx~a ~ 0(1). This in 
turn implies that the scattering is strong for projectiles with transverse sizes r > 1/Qs- For 
instance, the dipole scattering amplitude {Tx^x2)y becomes of order one when \xi — X2\ ^ l/Qs- 
Then Eq. (3.19) implies that, for generic configurations at least, the quadrupole scattering 
becomes strong when at least one (which necessarily means at least three) of the six transverse 
distances rij = \xi — Xj\ is of order 1/Qs, or larger. Similar considerations apply to the higher- 
point correlations. In this regime, the Wilson lines cannot be expanded out anymore. Rather, 
they resum multiple scattering to all orders in the eikonal approximation. 

To correctly describe the high-energy evolution in the presence of gluon saturation and mul- 
tiple scattering, it is of course essential to keep the non-linear terms in the Balitsky-JIMWLK 
equations, so like {SS)y in the equation (2.16) for the dipole S'-matrix and {SQ)y in the r.h.s. 
of Eq. (2.17) for the quadrupole. In fact, these are precisely the terms responsible for the 
approach towards saturation in the gluon distribution and towards unitarity in the scattering 
of the projectile. Accordingly, in the transition regime towards saturation/unitarity {i.e. for 
k± ~ Qs(Y)), one has to deal with the whole, infinite, hierarchy of coupled evolution equations: 
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no simple mean-field approximation (like a Gaussian) is possible in that regime. However, the 
situation drastically simplifies deeply at saturation {k± <C QsiX)), where the only role of the 
non-linear terms in the equation is to forbid further evolution — or, more correctly, to limit the 
transverse phase-space for the high-energy evolution: gluons with soft momenta k± <^ Qs{Y) 
can (almost) not be emitted anymore, meaning that domains separated by transverse distances 
r ^ l/Qs(Y) evolve independently from each other. This leads to considerable simplifications 
in the Balitsky-JIMWLK equations, which can be most directly recognized by inspection of the 
projectile evolution. 

For multi-partonic projectiles which are such that all the interparticle separations rij are 
much larger than l/Qs{Y), the associated ^-matrices are very small (close to zero) — the more 
so the larger the number of partons. Roughly speaking, and up to subtleties related to the l/N^ 
corrections to which we shall later return, a 2-dipole projectile scatters more strongly than a 
single-dipole one, {SS)y ^ {S)y, a projectile made with a dipole plus a quadrupole scatters 
more strongly than the quadrupole alone, {SQ)y «C {Q)y, etc. When this happens, the 'virtual' 
terms dominate the evolution, whereas the 'real' terms can be simply replaced with a lower 
cutoff ~ l/Qs{Y) on the transverse separation \z — Xi\ between the newly emitted gluon at z 
and any of the preexisting partons at Xi. Once this is done, the resulting evolution equations 
are linear and hence admit a Gaussian solution. This is of course related to our previous 
observation in Sect. 2.3 that the only effect of the 'non-linear terms' (Wilson lines) within -ffvirt 
is to transform 'left' color precessions into 'right' ones and thus ensure the symmetric expansion 
of the target field distribution in x~ . This also shows that, in this high density regime, where 
the Wilson lines cannot be expanded anymore and 'left' and 'right' functional derivatives have 
different mathematical consequences, it is essential to keep trace of the 'mirror' symmetry of the 
evolution, by using a symmetric Gaussian, as shown in (3.1). 

To render these considerations more precise and construct the corresponding Gaussian ap- 
proximation, we shall develop our mathematical arguments in two steps: (i) at large Nc, and 
(ii) at finite Nc- 

(i) Large Nc : Within the context of the large- A'c approximation, the prominence of the 
'virtual' terms in the approach towards the black disk limit is quite obvious and has been 
pointed out at several places in the literature [12, 39, 56, 57]. Specifically, the 'real' terms which 
survive at large Nc involve double-trace operators, which can be factorized to the accuracy of 
interest: {SS)y — {S)y{S)y, {SQ)y — {S)y{Q)y, etc. Then we can write e.g. 



Now, in equations like (2.16) or (2.17), the transverse position z of the emitted gluon is integrated 
over, so it can become close to one of the external points Xi, in which case Eq. (3.27) does not hold 
anymore. However, in the high density regime under consideration, such special configurations 
are disfavoured by the phase-space for the transverse integration. Namely, assuming \xi — Xj\ ^ 
l/Qs{Y) for all the pairs (i, j), one can check that the integrals over z receive their dominant 
contributions from points relatively far apart from all the external points, which satisfy 




when 



z-x^\ » 1/QsiY). 



(3.27) 



l/<5s < \Z - Xi\ < \Xi-Xj 



(3.28) 
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Indeed, the contribution of such a range is enhanced by the large transverse logarithm 
1 /" dz"^ 



„ , M^^^^z ^ / " ^ = ln[{x,-x,fQl] . (3.29) 



Hence, to leading logarithmic accuracy in the sense of Eq. (3.29), one can indeed neglect the 
'real' terms in the Balitsky-JIMWLK equations at large Nc, as anticipated. 

(ii) Finite Nc : The physical argument at finite Nc is the same as at large Nc except that, 
now, one has to take into account the fact that the evolution described by the 'real' terms truly 
corresponds to the emission of a gluon, and not just to the splitting of, say, one dipole into 
two dipoles. When this new gluon is sufficiently soft, in the sense that \z — Xi\ > l/Qs{Y) for 
any i, its emission leads to a partonic system with a wider distribution of color charge in the 
transverse plane, which therefore interacts stronger with the target than the original projectile. 
But in order to rigorously justify this, one needs to actually estimate the S-matrix for, say, a 
quark-antiquark-gluon (qqg) system deeply at saturation and show that this is indeed much 
smaller than the iS-matrix of the dipole (qq). To appreciate how subtle this is, let us recall that, 
when rewriting the 'real' terms in terms of Wilson lines in the fundamental representation (as 
customary in the Balitsky-JIMWLK equations), one generates single-trace pieces proportional 
to l/Nc, which by themselves count on the same footing as the 'virtual' terms near the unitarity 
limit. For instance, the contribution of the 'real' terms to the dipole equation involves the 
following expectation value (cf. the second line in Eq. (2.13)) 

SxizS *^^i^2^^ ' (3.30) 

where one may naively think that the second, single-trace, term dominates over the first one 
when all the transverse separations are much larger than l/Qs(Y). As another example, we 
show here some 'real' terms from the evolution equation (2.17) for the quadrupole, namely 
those arising when acting with .ffreal on the two quarks at a^i and x-^ : 

{H,..iQx^„,)y = -^:^J^Mx,x,z{vf[tT{^^^^^ 



/ ■^^xiXjz\^Szx2Qxizx3X4 ~\~ Szx4QxiX2X:jZ ^2 QxiX2X:jX4 j ^ ■ 

(3.31) 



' c 



where the second line follows from the first one after using the Fierz identity (2.15). Once again, 
one may think that the last term in Eq. (3.31), proportional to {1/N'^){Q)y , is the dominant 
term for large transverse separations l/QgiY) (and for finite Nc). If that was indeed the 
case, there would be a mixing between 'real' and 'virtual' terms deeply at saturation, which 
would prevent a Gaussian approximation (since the latter could not accommodate the 'real' 
terms beyond the BFKL approximation). 

The situation becomes even more confusing if one recalls that, in the equations obeyed by 
the single-trace observables, the terms subleading at large Nc precisely cancel between 'real' 
and 'virtual' contributions. In view of this, one may be tempted to argue that the finite-A'^c 



- 24 - 



corrections are totally irrelevant. But that would be wrong, since there is no similar cancelation 
in the equations obeyed by the multi-trace operators, like {SS)y or {SQ)y- 

What 'saves' the Gaussian approximation, is the fact that, in spite of appearance, the 
single-trace components in equations like (3.30) or (3.31) do not dominate over the respective 
double-trace ones, but merely subtract fake 'single-trace contributions' from the latter, that 
have been artificially introduced via the Fierz identity. That is, the expression in the first line 
of Eq. (2.13), which involves an adjoint Wilson line and describes a qqg system, vanishes very 
fast in the approach towards the black disk limit, where it is suppressed with respect to the 
corresponding 'virtual' term {Sxix2)y- But this is not the case for the 2-dipole S'-matrix in 
the second line of Eq. (2.13), which in that regime approaches to {1/N^){Sxix2)y- A similar 
discussion refers to Eq. (3.31): deeply at saturation, the observable in the first line, which 
describes a qqqqg partonic system, is suppressed compared to the respective 'virtual' terms, 
that is, the quadrupole and the pair of dipoles. 

In order to demonstrate this while dealing with an infinite hierarchy, we shall provide a 
self-consistent argument. That is, we start by assuming that the JIMWLK evolution deeply at 
saturation is controlled by -ffvirt alone and we prove that, under this assumption, the 'real' terms 
in Eq. (3.30) and Eq. (3.31) vanish exponentially faster than the respective 'virtual' terms in the 
vicinity of the unitarity limit. We shall give the details of the proof for the dipole evolution, i.e. 
for the operator in Eq. (3.30), and then briefly discuss its generalization to the quadrupole and 
higher n-point functions. In this context, by '-ffvirt' we mean, of course, the first two terms in 
the JIMWLK Hamiltonian (2.2) together with the phase-space restriction \z — Xi\ ^ l/Qs{Y) 
as introduced by the 'real' terms. That is, we work in the leading-logarithmic approximation 
in Eqs. (3.28)^(3.29), which enables us to write 



' virt 



-^£l„[(„-«)^Q2(y)] + (3.32) 



to the accuracy of interest. 

So, let us calculate the action of -ffvirt on the combination of the operators appearing in 
Eq. (3.30). This action on the second term has been already computed in Eq. (2.12), that we 
here rewrite for convenience as 

H^Yl-l Sx]^X2 ~ 7, aTo ( ^ Art ) / ■^X\X2wSxiX2-I (3.33) 



c 



with the integral over w understood in the sense of Eq. (3.29). Now, when both derivatives act 
on the same (either the first or the second) dipole of the first term in Eq. (3.30), we get the 
following, 'diagonal', contribution 

H^/ixt SxizSzx2\f^l^g ~ ^1 iV?^) j ^"'^"'i^'"' ~^ ^^X2Zva)SxizSzx2^ (3.34) 

and when they act on different dipoles we find the cross term 

-f^virt »S'a;jz5'2;a;2 IctoSS 27J-_/\r2 J i-^^XlZW ~^ ■^^X2Z'W ■^^XlX2'w){SxizSzX2 SxiX2^- (3.35) 



- 25 - 



Putting everything together we arrive at 



SxizSz 



+ M X1X2W 

27r ,„ V 



SxizSzX2 iy2 SxiX2 

c 



(3.36) 



where it is crucial to notice that the operator of interest has been reconstructed in the r.h.s. of 
the equation. It should be clear from the above derivation that this would have not happened 
without the subtraction of the l/A^^-suppressed dipole. By assumption, the above equation 
describes the approach towards unitarity of the 'real' piece in the evolution equation (2.16) for 
the dipole S'-matrix. This should be compared to Eq. (2.12), which describes the corresponding 
approach for the 'virtual' piece (the dipole itself). Clearly, the kernel in Eq. (3.36) is 'twice 
as large' than that in Eq. (2.12), showing that, deeply at saturation, the expectation value of 
the 'real' operator in Eq. (3.30) vanish exponentially faster than the 'virtual' term oc {Sx^x2)y- 
Hence, the latter dominates in the evolution equation and in this regime, as anticipated. 

This self-consistent argument can be generalized to higher-point correlators, as we now 
show for the operator 



SxizQ 



ZX2XgX4 



Qxia 



(3.37) 



which appears in Eq. (3.31) and counts for the evolution of the quadrupole. Acting with -ffvirti 
we see that the only new element appearing, when comparing to Eq. (3.36), is operator mixing. 
Indeed, one finds that we also need to consider the operators 



and 



Q Q c* c* c* 

^XlZ)^ ZX2^X:iXi JLJ2 ^XlX2^XzX4^1 

c 



(3.38) 



0(6) _ Q C 

'JzX2X\ZXgX4 'JxlX2^X-gX4,^ 



(3.39) 



plus permutations of all the operators appearing in Eqs. (3.37), (3.38) and (3.39). Without going 
into too much detail, one understands that the action of i/virt on the above operators leads to 



virt 



SxizQ ZX2XZX4 Qx\X2XZXi 

5^ Q Q C* C* 

xizi~'zx2'~'x:iX4 ^2 '~'xiX2'~'x[iX4 



s 



(6) 



ZX2XlZXgX4 



SxiX2 ^x:i 



3x4 



SxizQ ZX2X3X4 QxiX2X3X4 

Q Q Q ^ Q Q 

'JxlZI~'zX2'~'xgX4 '~'xiX2'~'X[iX4 



:(6) 



■'ZX2XIZXSX4 



SxiX2 



X3X4 



(3.40) 



where the elements of the 3x3 matrix M are proportional to —ajlTx times an integral over 
w of linear combinations of the dipole kernel. The counting is such that the integrand in the 
diagonal elements is the sum of three dipole kernels which enter all with a plus sign, plus terms 
proportional to l/A'^^ (in analogy with Eq. (3.36)). Furthermore, the integrand in the non- 
diagonal elements is the sum of dipole kernels with equal number of plus and minus signs, plus 
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again terms proportional to Clearly, the diagonal components are those which control 

the approach towards the black disk limit and they are larger than those which control the 
corresponding evolution for the 'virtual' terms in Eq. (2.17), that is {Q)y and {SS)y- 

Incidentally, the above argument also shows that the two operators in Eqs. (3.38) and 
(3.39) vanish faster than the quadrupole and the 2-dipole system in the approach towards 
unitarity. This is interesting since these are precisely the 'real' terms in the evolution equation 
for {Sx^x2Sx3X4)y , whereas (Q)y and {SS)y are the corresponding 'virtual' terms. So, we 
have also demonstrated the property of interest (the dominance of the 'virtual' terms deeply 
at saturation) for the evolution of a system of two dipoles with arbitrary coordinates. We are 
confident that a similar proof applies to the higher-point (single-trace or multi-trace) correlation 
functions. 

It is furthermore instructive to check these arguments via explicit calculations within the 
Gaussian approximation (3.1). Via methods to be described later, this yields e.g. [37] 



1 \ iVc 1 

'S'a;i£C3'S'a;3a;2 ~ 3x^x2 J ^ = 



(5'a;ja33)y (»S'a;3a;2 )y 
{8x1x2) Y 



{Sxxx3)y{Sx3X2)y-, (3-41) 



where we have assumed that {Sx^xj) and {8x1x38x3x2) are equal to 1 as an initial condition, to 
simplify writing. This formula makes it clear that the operator in the l.h.s. vanishes, roughly, 
as a 'dipole squared' in the approach towards the unitarity limit. A corresponding argument for 
the operator (3.37) which enters the evolution of the quadrupole will be given in Sect. 4.4. 

We thus conclude that the JIMWLK evolution deeply at saturation is indeed correctly 
described by the 'virtual' Hamiltonian in Eq. (3.32). When acting on operators built with 
Wilson lines, the two terms in i^virt amount to 'left' and 'right' Lie derivatives, in the sense of 
Eq. (2.25). So, clearly, the Hamiltonian (3.32) is of the 'symmetric Gaussian' form in Eq. (3.4), 
with the following kernel 

jy{u, v) = ^In [{u - vfQl] =^ jY{k) = ^ . (3.42) 

This applies for k± <C Qs{Y) and is recognized as the 2-dimensional Coulomb propagator. In 
turn this implies that the charge-charge correlator Ay(fc) = A;^7y(fc) vanishes like k'j_ when 
fc_L — 7- 0, which is the expression of color shielding due to gluon saturation [38, 58]: the average 
color charge squared vanishes when integrated over a transverse area ^ 1/Q'^{Y). 

Notice that in some previous versions of the mean field approximation [36, 38, 39], one 
has assumed that the JIMWLK Hamiltonian takes an even simpler form in the vicinity of the 
black disk limit, namely it reduces to the first term in Eq. (3.32), which involves the 'left' 
derivatives alone. That simplification was motivated [39] by a 'random phase approximation', 
which assumed that, in the strong field regime deeply at saturation, all the Wilson lines within 
the Hamiltonian are rapidly oscillating and thus average out to zero. As shown by our present 
manipulations, this argument is qualitatively correct, but only for the 'real' terms (the last 2 
terms) in the JIMWLK Hamiltonian. 

To summarize the arguments in this section, the JIMWLK evolution in the two limiting 
regimes — the weak-scattering regime at low gluon density and the approach towards the 
black-disk limit deeply at saturation — can be properly encoded, for any value of Nc, into a 
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symmetric Gaussian weight function of the type (3.1). In turn, this Gaussian is tantamount to 
the functional evolution equation shown in Eq. (3.4), or to the following, mean field, Hamiltonian: 

HMFA = -\j iY{u,v)[l + VlV.y'^^ = ^^ f 7y(«,^)(JE.J£„ + J]^.JU- (3-43) 

This has the same operator structure as the 'virtual' piece of the JIMWLK Hamiltonian, but 
with a different, y-dependent, kernel, which is essentially the 2-point function of the target 
color field. This kernel interpolates between the solution to the BFKL equation (3.26) at small 
transverse separations |m — i>| <^ l/QgiY) and the Coulomb propagator (3.42) at relatively 
large distances |u — i»| ^ \/Qs(Y). Remarkably, the kernel is independent of Nc, in agreement 
with the corresponding property of the 'dipole' kernel in the JIMWLK Hamiltonian. (This can 
be checked e.g. on Eq. (3.42) and on the relation (3.24) between this kernel and the dipole 
amplitude at weak coupling.) Any smooth function 7y(it, v) with the correct limiting behaviors 
can in principle be used to define the Gaussian; indeed, different such functions can differ from 
each other only in the transition region around Qs , which is not under control within the present 
approximation. In practice, however, a proper choice for the kernel is probably essential in order 
to achieve a good global accuracy. In the following section, we shall propose two such choices. 



4 Evolution equations in the Gaussian approximation 

In this section we shall describe two methods for constructing smooth, global, expressions for 
the kernel 7^(^,1?), which are equivalent with each other to the accuracy of interest. Then 
we shall derive the evolution equations associated with the mean field Hamiltonian (3.43), first 
for generic Nc (in Sect. 4.2), then at large Nc (in Sect. 4.3). As before, we shall mostly focus 
on the evolution of the dipole and of the quadrupole. In the large-iVc limit we shall recover 
the equations previously proposed in Ref. [12]. At finite Nc, the general equations are more 
complicated, but explicit solutions will be presented for special configurations in Sect. 4.4. 

4.1 Self consistent constructions of the kernel 

Within the Gaussian approximation, it is always possible to trade the kernel 7y(a;i, X2) for the 
dipole ^-matrix {Sx-^x2)y> for which it is easier to construct global, smooth, approximations 
in practice. The expression of {8x1x2) Y in the Gaussian approximation is well known in the 
literature and will be rederived, for completeness, in the next subsection. Here it is preferable 
to work with the corresponding evolution equation, which reads 

^^%a>i: = -2g^Cnly{a^i,x,){S^,x2)y (4-1) 

for a dipole in an arbitrary representation R of the color group. Hence, if one disposes of a 
numerical solution to the JIMWLK equation, like in Refs. [11, 59, 60], then one can use the 
respective estimate for the dipole S-matrix, say, in the fundamental representation, together 
with Eq. (4.1) to deduce a corresponding estimate for the kernel. 

In practice, solving the full JIMWLK evolution is quite tedious, so it is customary to rely on 
its large-iVc approximation (for the dipole evolution), namely the BK equation. This is equally 
good for the present purposes (including at finite Nc) since, as noticed at the end of the previous 
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section, the kernel 7y(a;i,a;2) is independent of Nc- Hence, its limiting behaviors are correctly 
reproduced by the large- A'c version of Eq. (4.1). The latter implies 



ain(5^?^l\^ 

g'N, ^Yix,,X2) = ^-^^ , (4.2) 

with {S^^2)y denoting the solution to the BK equation with an initial condition itself evaluated 
at large Nc- (This is important in order to preserve finite- A'^c accuracy in the limiting kinematical 
domains where Eq. (4.2) is correct as it stands for any value of Nc-) For instance, within the 
MV model, this is provided by Eqs. (3.11)-(3.12) with Cp ^ iVc/2. The function {S^^X2)y 
can be obtained either as an exact, numerical, solution to the BK equation, or as an analytic 
approximation to it with the correct limiting behaviors. 

Although correct to the accuracy of interest, the construction in Eq. (4.2) may look a bit 
unesthetic at a conceptual level, as it requires an input — the solution of the BK equation — 
which seems external to the Gaussian approximation. However, as we now explain, Eq. (4.2) 
can be also understood as a self-consistency condition internal to the mean field approximation 
[37]. Specifically, let us start with the first equation in the Balitsky-JIMWLK hierarchy, that is 
Eq. (2.16) for the dipole, and evaluate all the expectation values there with the Gaussian weight 
function (3.1), that is, by using Eq. (4.1) with Cr = Cp together with Eq. (3.41). This leads to 
an equation for the kernel 7y(a;i,a;2) which is precisely equivalent to solving the BK equation 
and then computing the kernel according to Eq. (4.2) (see Ref. [37] for details). 

This procedure, which in Ref. [37] has been dubbed 'the Gaussian truncation', is of course 
not unique: one can similarly start with any equation in the Balitsky-JIMWLK hierarchy, 
compute all the expectation values there with the Gaussian weight function (3.1), and thereby 
transform the original equation into an equation for 7y(iCi,a;2)- A different self-consistency 
condition, which in practice is not more difficult to use than Eq. (4.2), has been originally 
proposed in Ref. [36]. It amounts to requiring the mean-field Hamiltonian (3.43) to coincide 
with the JIMWLK Hamiltonian (2.2) on the average : 

^ ^ M^..{1 + KlK - Vl% - V^X)t = \ 7y ^) (1 + (4-3) 

The average here refers, of course, to the Gaussian weight function, which implies that both the 
l.h.s. and the r.h.s. in the above equation are proportional to b""^ . Introducing the ^-matrix for 
the gluonic dipole operator 



stx2 ^ 'iv(yi^Kj = [nI Sx,x2S.2x, - 1) (4.4) 

which it related to the respective fermionic operator as shown in the second equality above, and 
multiplying Eq. (4.3) by S"'^, we can rewrite the latter as 

= ^ ^ A. V / Mnvz{l + St - St - St)Y- (4-5) 



1 „ ~ . 1 



uv I Y 



This relation immediately implies 7y(it, u) = because of the corresponding property of the 
dipole kernel M.uvz- It furthermore implies that 7y('U, t;) is symmetric under ti -f-^- v, because 
so is the gluonic dipole S'-matrix, as obvious from its definition (4.4). 
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The self-consistency condition (4.5) is most conveniently written as an equation for (5*^)^: 
by using Eq. (4.1) in the adjoint representation {Cr = Ca = -/Vc), one finds 

QY ~ TT '3313322 cA \ \ ^ '^ZX2 "^xiXi ^/y V^-'^) 

The initial condition should be taken from the MV model applied to a gluonic dipole, that is, 
Eqs. (3.11)-(3.12) with Cf^Ca = N^. 

To summarize, by using either Eq. (4.2) together with the solution to the BK equation, or 
Eq. (4.1) with R = A together with the solution to Eq. (4.6), one obtains two global approxi- 
mations for the kernel 7y (a^i, 0^2), which can differ from each other only in the transition region 
around saturation. In principle, these two approximations are equivalent with each other to the 
accuracy of interest. In practice, any (numerical) difference between them will have an impact 
on the calculation of the higher n-point functions, to be described in the remaining part of this 
section. It is likely that one can optimize the reliability of this whole scheme by choosing a 'good' 
approximation for the dipole ^-matrix used to compute the kernel, like the exact, numerical, 
solution to the JIMWLK equation, or to the BK equation at least. 

4.2 Mean field equations for the dipole and the quadrupole 

The evolution equations associated with the mean-field Hamiltonian (3.43) are straightforward 
to obtain, by using the same techniques as for the JIMWLK Hamiltonian (2.2). Alternatively, 
given that i^MFA has the same operator structure as -ffvirti the mean-field equations can be 
directly inferred from the corresponding Balitsky-JIMWLK equations, by keeping only the 'vir- 
tual' terms in the latter and replacing everywhere the dipole kernel according to 
1 

8^ 

In doing that, one should be careful to restore all the 'virtual' terms in the Balitsky-JIMWLK 
equations, including those which may have canceled against similar 'real' contributions and hence 
are not manifest in the final equations (cf. the discussion in Sect. 2.2). Clearly, the resulting 
equations will inherit the relatively simple structure characteristic of the 'virtual' terms. They 
form closed systems of equations, which connect only correlation functions with the same number 
of Wilson lines (or external points) and are local in the transverse coordinates, meaning that 
they do not mix different transverse configurations. Note also that, although linear, these new 
equations respect unitarity by construction: the tame of the BFKL growth by the non-linear 
physics of gluon saturation is already encoded in the kernel of the Gaussian. The fact that the 
S-matrices for the various projectiles approach the right limit at strong scattering is ensured by 
the unitarity of the Wilson lines which appear within the respective operators. 

Consider first the dipole equation. Using the substitution rule (4.7) and the respective 
'virtual' term in Eq. (2.12), one finds 

d{S^^X2)Y ^ -2g'CFjy{x^,x,){S^,^,)Y . (4.8) 
This is easily solved to give 

fY 

(S.i=.2)y = e-r^("i'"^), TY{x^,X2) = 2g^CF dy^yix^^x^) , (4.9) 

J —00 



Muvz lYiu,v) . (4.7) 



z 
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where it is understood that Tyq is given by the MV model, cf. Eqs. (3.11)-(3.12). The corre- 
sponding expressions for a color dipole {S^^^2)y arbitrary representation R are obtained 
by replacing Cp — ?• Cr in the equations above (cf. Eq. (4.1)). In particular, in the 'BK- 
representation' in which the Gaussian kernel is computed according to Eq. (4.2), the dipole 
S-matrix in the Gaussian approximation and in a generic representation R of the color group 
is related to the solution {S^^^^)y to the BK equation via 



ln(5«.,>^ = ^ln(5^f.,>^. (4.10) 



From the discussion in Sect. 3, one expects Eq. (4.9) to have the correct limits at both weak 
and strong scattering, and it is instructive to explicitly check that. At weak scattering, one has 
{S)y = 1 — {T)y with {T)y ^ 1, so in particular one can replace {S)y ~ 1 in the r.h.s. of 
Eq. (4.8). Then the latter reduces to Eq. (3.24), with 7y(a?i,a?2) determined by the solution to 
the BFKL equation (3.26). This is indeed the expected result. At strong scattering, the kernel 
takes the form of the Coulomb propagator (3.42), in agreement with the limit |w — i;! » l/QgiY) 
of Eq. (4.7) (recall Eq. (3.29)). So, in this regime, Eq. (4.8) is identical to the 'virtual' part of 
the respective Balitsky-JIMWLK equation (2.16), which is the part that controls the approach 
towards the black disk limit. Let us study this approach in more detail. By using Eq. (3.42) 
within Eq. (4.9), one obtains 



27r^ lYJr-\ ./l/r.2 KI JvJrA 



Y - Ur)y = — ^-^ In^ {r'Ql{Y)) (4.11) 



where Ys{r) is the rapidity at which saturation is reached over a transverse size r (that is, 
Qs{Ys{t)) = 1/?"), is the logarithmic derivative of the saturation momentum, and we used 
In (r^Q^(y)) = U}6i{y — Ys{r)) for y > Ys{r). Eq. (4.11) holds in the leading logarithmic approx- 
imation in the sense of Eq. (3.29). The large- A'^c version of this result has already appeared in 
the literature [56, 57]. 

We now turn to the quadrupole. The corresponding evolution equation in the MFA is 
obtained as 

d{QxiX2X3Xi)Y = _ g'i(JpY^Y{xi,X2) +-fY{xz,X2) + -iY{xz,Xi) + ^Y{xi,XA)]{QxiX2X-iX4)Y 

[27y(a;i,a;3) + 27^(2:2, 2:4) - ^y{xi,X2) - 7y (2:3, 2:2) - 1y{xz 

,2^4) ^y{x1iX4)\{QxiX2X3,X4)y 



2N, 



- ^-^^[1y{xi,X2) +'^Y{X'i,Xi) -^y{xi,X^) - '^Y{x2,X4)]{Sxix^Sx3Xi)Y 

- ^-^[lY{xi,Xi) +-iY{x^,X2) - 7y(a;i,a;3) - ^Y{x2,X4)]{Sx^XiSx:,x2)Y ■ (4.12) 

Once again the BFKL limit is easy to check: when all the separations | much smaller 

than 1/QsiY), one can replace {Q)y ~ 1 and {SS)y ~ 1 in the r.h.s. of the above equation, 
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which then reduces to 
_d_ 
dY 



{QxiX2XiX4)Y ^ -2g Cf[7y{x2,X3) +Jy{xi,X4) +-fY{xi,X2) +7^(033, 034) 



7y(a;i,a;3) - -fY{x2,X4)] 



TxiX4 ~1~ Txj^X2 ^ Tx2 



x:jX4 



'^X2X4)Y 



(4.13) 



The second hne, which fohows from the first one after using Eq. (3.24), is the expected result 
for {Q)y at weak scattering, cf. (3.19). 

Let us also notice that, within the Gaussian approximation, the dipole and quadrupole S- 
matrices are invariant under charge conjugation, that is, under the exchange of the quarks with 
the antiquarks. More precisely, from Eqs. (4.8) and (4.12), and using the fact that the kernel 
'yY{xi,Xj) is symmetric, we easily deduce that 

{SxiX2)y — {Sx2Xi)y and {QxiX2XsX4}y — {Qx2XgX4Xi)Y ^ ('^•-^'^) 

so long as the above conditions already hold at Yq (as is the case within the MV model). 
Conversely, C-odd scattering amplitudes, like the odderon, cannot be accounted by the Gaussian 
approximation, as they would require non-Gaussian effects already at Yq (cf. the discussion at 
the end of Sect. 3.3). 

Returning to the full equation (4.12), we notice that this is consistent with mirror symmetry, 
as it should. (E.g. the r.h.s. is symmetric under the exchange xi o x^.) That would not have 
been the case^, had we considered a Gaussian Hamiltonian with only left derivatives, as proposed 
in Refs. [36-38]. In fact, without any further assumption, this equation implies that the evolution 
for the antisymmetric part of the quadrupole defined in Eq. (2.20) is closed and homogeneous. 
In turn, this means that (Q^**^™)y = in the MFA provided this condition was satisfied at Yq. It 
is furthermore clear that the quadrupole couples to the product of two dipoles, whose evolution 
in the MFA is in turn determined by 



d{SxiX2 Sx'^X4)y 



dY 

^2 



2g'^CF[7Yixi,X2) + Jy{x3, X4)]{SxjX2Sx3X4)y 



■[7y(cCi,£C3) + Jy{x2,X4,) - Jy{xi,X4) - -fY {x^, X2)]{Sx^x2SxaX4)Y 



2N, 



\-lY{xi,Xi) +^y{x3,X2) - ^y{xi,X3) - 7y(£C2,£C4)]((5a,ia 



{X1X4X3X2/Y ■ 



(4.15) 



Since the equations above involve {SxiX4Sx3X2)y and {QxiX4X3X2)y we need also to consider 
them with the their indices X2 and x^ interchanged. (Actually, {QxiX4X3X2)y coincides with 
{QxiX2X3X4 if one assumes mirror symmetry at Yq, but here we prefer to keep the discussion 
general.) Thus we arrive at a homogeneous system of first order differential equations 



_d_ 
W 



{QxiX2XsX4)y 
{QxiX4XsX2)y 
IX2^X3X4)y 

{Sx 



MY{Xi 



{Q X1X2XSX4) Y 
{Q X1X4XSX2) Y 

{SxiX2SxsX4)y 
{SxiX4Sx3X2)y 



(4.16) 



^To see this, let us assume for the sake of the argument the large-A^c hmit where the second term is absent. 
Then we find that in the first term only 'yvixs, X2) and 7y (cci, 0:4) are present and that the fourth term is absent. 
Clearly the aforementioned symmetry of the evolution equation is lost. 
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with My a 4x4 matrix. Its elements are proportional to (7^7^ (a;,, aij) accompanied by color 
factors and can be read from Eqs. (4.12) and (4.15). The difficulty that appears now is that 
one cannot solve Eq. (4.16) for a generic dependence of ^Y{xi,Xj) on Y, since in general the 
matrices My at different rapidities Y do not commute with each other, that is [My^jMyj] 7^ 0. 
(More precisely, one could write down a formal solution which involves the rapidity-ordered 
exponential of the mixing matrix, but we do not find that very useful in practice.) 

There are special cases where the rapidity integration in the equation above can be explicitly 
performed, leading to a simpler expression. The large-A'^c limit to be discussed in the next 
subsection is one such a special case. Another one is when the kernel 7y(a;j, Xj^^ is El separable 
function of Y and the transverse coordinates, plus an arbitrary function of Y : 

-fvixi, Xj) = hi{Y) g{xi, Xj) + h2{Y) . (4.17) 

This property is manifestly satisfied within the (homogeneous) MV model, as noticed after 
Eq. (3.9), and it is also approximately satisfied by the solution to the BK equation, at least in 
particular kinematical regimes: in the window for extended geometric scaling, where 7y(r) oc 
{r'^Ql{Y))"<'' with 7^ w 0.63 [49-51], and also deeply at saturation where 7y(r) oc ln[r^Qg(y)], cf. 
Eq. (4.11). Presumably, this is a reasonable approximation for all the dipole sizes. Furthermore, 
this is also fulfilled in some widely used dipole models, like the GBW model [61, 62], where 
ry(^) = f'^Qliy) The role of separability in simplifying the results of the high-energy 
evolution in a similar context has been previously recognized in Ref. [35]. 

Within such a simplified scenario, the y-dependence factorizes out from the mixing matrix 
and the resulting, y-independent, matrix M can be diagonalized. Then one can explicitly solve 
the system of equations. Within the context of the MV model this was done in Refs. [8, 31, 35]. 
In that case, one had to deal with a 2 x 2 system only, because the Wilson lines were constructed 
via 'left' iterations alone, cf. Eqs. (3.16)-(3.17). That is, the associated, effective, evolution 
equations were not explicitly 'mirror-symmetric', unlike the above equations (4.12) and (4.15). 
However, this symmetry is recovered in the final results, because the color charge distribution 
in the MV model in symmetric in x~ (cf. the discussion after Eq. (3.17)). Moreover, for a 
separable kernel, the final results for all the n-point functions depend only upon the integral 
/ ^Uly (a property which looks natural in the case of the dipole, cf. Eq. (4.9), but which in 
general requires separability). Hence, the results obtained in [8, 31, 35] within the MV model 
can be transposed to a more general Gaussian which is still 'separable', by simply replacing 
the kernel in the final formulae according to Eq. (4.17). We shall not write here the respective 
general results, but refer to [8] for {Qx^x^xiXi)Y and to [35] for {SxiX2Sxsxi)Y ■ 

But separability is not always needed in order to obtain explicit solutions at finite N^- for 
special configurations of the 4 external points Xi in the transverse plane, the matrix My in 
Eq. (4.16) may happen to simplify independently of the structure of the kernel ^yixi^Xj). As 
an example consider the evolution of the 2-dipole 5-matrix {Sx-i^xzSx3X2)Y , in which the quark 
in one dipole and the antiquark in the other dipole are located at the same point 2:3. The 
expression of this correlator in the Gaussian approximation has been shown in Eq. (3.41) and 
we would like to check that here. By identifying X2 with x^ in Eq. (4.15) and then relabeling 
X4^ as X2 for convenience, we see that the two quadrupole operators in the r.h.s. there reduce to 
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single dipoles and then this equation decouples from the evolution of the quadrupole 



d{Sx 

dY 



+— [7y(a;i,a^3) + 7y(a;3,a;2) 



■7y (031,032)] } 

■7Y(a;i,a;2)] {SxiX2)y- 



(4.18) 



The last, inhomogeneous, term in the r.h.s. is particularly interesting, as it describes a process 
in which the two dipoles Sxixg and Sx-jX2 having one common leg merge with each other into 
a single dipole SxiX2- This process is suppressed at large Nc (as expected, since dipoles cannot 
merge with each other in the large Nc limit [27]) but for finite Nc it controls the approach 
towards the unitarity limit, since a single dipole scatter less than a system of two dipoles. 

Using Eq. (4.8) to express ^Y{xi^Xj) in terms of the logarithmic derivative of {SxiXj)Y we 
can rewrite the above equation as 



d{Sxix^ Sx'^X2}y 

dY 



9 {SxiX3)y {Sx3X2)y 

dY ^ 



l+e 



{Sxi 



X2)y 



{Sx^x-iSx-j,X2)Y 



9 {SxiX3)y{Sx3X2)y 

dy {Sx,x2)y 



{SxiX2)Yi 



(4.19) 



where we temporarily defined e = l/{Nc — 1). This is a first order inhomogeneous linear 
differential equation which can be readily solved, with the result shown in Eq. (3.41). Some other 
special configurations, which in particular involve the quadrupole, will be studied in Sect. 4.4. 

The generalization of the above considerations to an arbitrary n-point function like Eq. (2.9) 
is straightforward. For instance, the mean-field version of the equation obeyed by the sextupole 
5-matrix (S^^)) Y can be inferred from the results in Appendix B of Ref. [12]. 

4.3 Quadrupole evolution at large Nc 

In the large- A'^c limit, the hierarchy generated by the Gaussian approximation drastically sim- 
plifies: it reduces to a triangular hierarchy, in which the equations can be successively decoupled 
from each other and explicitly solved. Let us illustrate that on the example of the 4-point 
functions — the quadrupole {QxiX2xax4)Y and the 2-dipole system {SxiX2SxiX4)Y — which in 
general mix under evolution, as shown in Eq. (4.16). At large N^ one can ignore the last two 
lines in the r.h.s. of Eq. (4.15), meaning that the 2-dipole system evolves independently of the 
quadrupole. The corresponding entries in the mixing matrix in Eq. (4.16) are now equal to zero, 
so that this matrix becomes triangular, as anticipated. Specifically, Eq. (4.15) reduces to 



d {Sx]^X2 ^X3X4)y 

dY 



-g^ Nc[-iY {xi,X2) + 7y(a;3, Xi)]{Sx^x2Sx-iXi)Y , 



(4.20) 



which is immediately solved as 



{SxiX2Sx3X4)y — 6 



ry, Yh (3:1 .3:2 ) -ry,y„ (tcs ,3:4 ) 



{SxxX2Sxs,X4)Yoi 



(4.21) 
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and where we have defined 



(4.22) 



As expected, this implies that the 2-dipole S-matrix factorizes at large Nc provided it did so in 
the initial conditions at Yq: 



{SxiX2SxsX4)y — {SxiX2)y {SxsX4')y ■ 



(4.23) 



Consider now the evolution of the quadrupole: by keeping in (4.12) only the leading terms 
at large Nc (which in particular means using the factorization property (4.23)), one finds 



_d_ 
W 



{QxiX2X^X4)y — TT' 

- [yY{xi,X2) + Jy{x3, X2) + Jy{x3, X^) + 7y(a;i, 0^4)] {QxiX2X3X4)y 

g'Nc 



g'Nc 



[yY{xi,X2) +Jy{x3,X4) - Jy{xi,X3) - 7y (032, 2:4)] {Sx-4X2)y {SxsX4)y 
[yY{x3,X2) +Jy{xi,X4) - Jy{xi,X3) - 7y (032, 2:4)] {Sx-sX2)y {SxiX4)y 

(4.24) 

which is an ordinary, first order, inhomogeneous differential equation. The dipole S-matrix is 
given by (4.9) with Cp ~ Nc/2 (as appropriate at large Nc) and acts as a source for the evolution 
of the quadrupole. As explained in Sect. 4.1, in practice it is preferable to view Eq. (4.9) as 
a definition of the Gaussian kernel in terms of the dipole (S-matrix, since the latter is a more 
directly relevant physical quantity. By using Eq. (4.2) to express ^Yixi,X2) in terms of the 
logarithmic derivative of {SxiX2)y we can rewrite Eq. (4.24) as 



d{QxiX2X't,X4)Y 1 



dY 



d_ 
dY 



111 ( 5'a3 a; 2 ) y ( 'S'tc 3 a; 2 ) y ( 'S'a; 3 a; 4 ) y ( 'S'a; 1 a; 4 ) y 



{Qxi 



X2XSX4/Y 



+ 



(5'£cia;3)y(<S'a,2a;4)y 

_d_ ('S'a;i3:4)y('S'a;3a;2)y 
(5'a;ia;3)y (5'a;2a;4)y 



( 'S'a; 1 a; 2 ) y ( "5*0; 3 £C 4 ) y 
('S'a;ia;4 ) y (5*033 332 ) Y ■ 



(4.25) 



The general solution of this equation is easily found as [12] 



(Qa;i3;2a;3a;4)y — y (»S'a;ia;2)y (5'a;33;2)y ('S'a;33;4)y (»S'3;ia;4)y 

xiX'i) y{Sx2X4) y 



{Q X1X2X3X4/Y0 



+ 



dy 



{SxiX2)Yo {SxsX2)Yo {Sx-sX4)Yo {SxiX4)Yo 

d {^xiX2)y {Sx ■iX4)y + {8x1X4) y{SxsX2)y 



Yo 



{Sx iX2)y {Sx ■sX2)y{SxsX4)y{SxiX4)y 



dy 



{Sxixs)y {Sx 2X4)y 



(4.26) 



This solution is already explicit, but it takes an even simpler form if one assumes separability, 
in the sense of the discussion after Eq. (4.16). In that case, the integral over y in Eq. (4.26) can 



- 35 - 



be exactly performed, to yield 

I /\ \ _ -^XlX2X3X4, \ /Q \ I '^XlX4XiX2 /Q \ /Q \ 

\WxiX2X[iX4/Y — J- \*-'a;ia32 /y \'-'a;8a;4/y ~r ^ \>->xzX2 /Y \>~>xiX4/Y 

^XiX2X4Xs ''-'X1X4X2X3 



+ 



{SxiX2)y {SxsX2)y {Sx3X4)y{SxiX4)y 

{Sx-iX2}yo ('S'a;3a;2 ) Vo {^X3X4 1X4/Y0 

i(\ \ ^'x\X2Xi,X4 iQ \ IQ \ ^x\X4XiX2 iQ \ I C \ 

\^x\X2XzX4lY(i J \>^xxX2lYa\>^xzX4lYQ j \>^xzX2lYa\^x\X4lYa 

^XiX2X4Xs ^X1X4X2X3 

(4.27) 

where we have denoted 

LxiX2X3X4 = '^Y,Yoixi,X2) + ^Y,Yoix3,X4:) - Ty^yp (aJi , 2:3) - Ty,Yo{^2, X^) ■ (4.28) 

Notice that the function L also depends upon Y and Yq, but because of separability the ratio 
between two L's is a function of the transverse coordinates alone. Eq. (4.27) depends upon the 
kernel 'yy{xi,Xj) only via its integral over y. This is a consequence of separability, as already 
noticed in Sect. 3.2 in the context of the MV model. As a matter of facts, Eq. (4.27) is quite 
similar to the corresponding expression in the MV model [3, 8] and it becomes formally identical 
to it once we assume an initial condition of the MV type. Specifically, Eq. (4.27) reduces to 

{QxiX2X3X4 )y — ^ ^ {SxiX2)y{Sx3X4)y + r ^ ^ ^ {Sx3X2)y{SxiX4)y, (4.29) 

-'^a;ia;23;43;3 '^xiX4X2X3 

provided this functional relation is already satisfied at Yq, as is indeed the case in the MV 
model and for large Nc [3, 8]. Note that there is an alternative way to deduce Eq. (4.29) from 
Eq. (4.27), which makes no reference to the MV model. Namely, if one assumes Eq. (4.27) to 
capture the whole evolution from 1" — )■ — oo (where the Wilson lines reduce to the unit matrix) 
up to the rapidity Y of interest, then one can use {Q)yo — ^ 1 and {S)yq — )• 1 for Yq — )• — oo to 
check that the expression within the square brackets in Eq. (4.27) vanishes for that particular 
initial condition. 

In general, i.e. without assuming separability, one expects the ratio of two L's to depend 
very weakly on Y. If so, it might be still a good approximation to use the simpler formula (4.29) 
for the quadrupole rather than the general one in Eq. (4.26), which is more involved. For that 
purpose, the function L in Eq. (4.29) should be defined by Eq. (4.28) with Yq — ?• — oo, and hence 

rY,Yo{xi,X2) rY{xi,X2) = - ln{Sx^x2)Y- (4.30) 

Given a smooth approximation for {S)y, such as the numerical solution to the BK equation, 
Eq. (4.26) (or (4.29)) provides a correspondingly smooth approximation for {Q)y, which is 
guaranteed to be correct whenever all the transverse separations r^j = \xi — Xj\ are either 
much smaller, or much larger, than l/Qs{Y), and for large Nc- On the other hand, the present 
approximations are strictly speaking not under control in the transition region around saturation 
{rij ~ l/Qs{Y)), nor for very asymmetric configurations, such that some distances r^j are much 
larger than 1/Qs while the other ones are much smaller. Some very asymmetric but relatively 
simple configurations will be discussed in the next subsection, directly for finite Nc- 



- 36 - 



4.4 Special configurations at finite Nc 

In this subsection, we shall study some special configurations of the 4-point function and the 
6-point function in the transverse plane, which because of their degree of symmetry allow for 
explicit, and relatively simple, solutions without additional assumptions like separability or 
large- A^'c. 

First we shall consider the class of configurations introduced in [12] for which the only 
constraints are ri3 = ri4 and = r24. For example, three such configurations are shown 
in Figs. 3. a, 3.b and 3.c. From these figures, it should be clear that there is a high degree 
of variability (concerning both shapes and sizes) within this particular class. By using the 
constraints aforementioned, it is straightforward to see that Eqs. (4.12) and (4.15) reduce to 

- ■^j;^bY{xi,Xs) +'yY{x2,Xi) - -fY{xi,X2) - 'yY{x3,X4)]{QxjX2X3X4,)Y 

[yY{xi,X2) +'yY{x3,Xi) --/y{xi,X3) - 'yY{x2,X4)]{Sx^x^Sx3X4)Y (4.31) 



2N, 



and 

d{S^,X2SxsX4)Y _ .2g'CF[^y^X,,X2)+-fY{x3,X2)]{S^,x2Sx3X4)Y. (4.32) 



dY 

Thus, for such configurations, the evolution of a system of 2 dipoles decouples from that of the 
quadrupole even without invoking separability. By also using Eq. (4.8) for the dipole S'-matrix 
in the Gaussian approximation at finite Nc, we can easily solve Eq. (4.32) to find 

/Q Q \ — {^xiX2Sx3X4)Yo / d \ / Q \ (A 
\JxiX2'->xaX4/Y — - - \'JxiX2/Y\'->X3X4/Y {^.OO) 

\!^xiX2)Yo \^xaX4)Yo 
which simplifies furthermore to 

{SxiX2^X[iX4)Y — {SxiX2)y {8x3x4) Y 1 (4.34) 

provided the latter holds in the initial condition at Yq (as is indeed the case within the MV model, 
as one can similarly check). Then it is clear that Eq. (4.31) can be solved as an inhomogeneous 
first order differential equation and it gives 

{QxiX2X3X4)y — r r {SxiX2)y {Sx3X4)y ~I~ {{QxiX2X3X4 )Yo - {Sx4 X2^X3X4)Yo] 

\JxiX2)Yo \>^X3X4)Yo 

2(nI~1) 



{Sx4X2)y {Sx3X2)y {Sx^X4)Y {Sx4X4)y I {SxiX4)y {Sx3X2)y {SxiX2)Yo{Sx3X4)Yo 



{SxiX2)Yo{Sx3X2)Yo {Sx3X4)Yo {SxiX4)Yo \ {SxiX2)y {Sx3X4)y {SxiX4)Yo{Sx3X2)Yo 



(4.35) 



Again, when the initial condition is given by the MV model (where Eq. (4.36) below holds 
indeed, as one can explicitly check), the above becomes very simple: 

{QxiX2X3X4)y — {SxiX2)y {Sx3X4)y • (4.36) 



-37- 



Xi 

• Xi 




(a) (b) (c) (d) 

Figure 3. (a), (b) and (c) correspond to special configurations of the quadrupole for which ri3 ~ ri4 
and r23 = r24. In (a) ri2 = and all Vij of the same order, in (b) r^i <C ri2 ^ rn and in (c) 
Ti2 ^ fsi ^ ''14 • The average values of Q and Sq depend only on the distances depicted by straight lines. 
Figure (d) corresponds to the line configuration for the operator Sq- 



This result is truly remarkable: within the class of configurations at hand, the quadrupole 
factorizes into two dipoles independently of how small or large the various distances r^j are, and 
for any Nc, so long as the two constraints ri3 = ri4 and = r24 are satisfied. It would be 
interesting to check this factorization via numerical solutions to the JIMWLK equation, as a 
non-trivial test of the present MFA. 

Now let us proceed to our second example and consider the operator 

Q.,.,.,.J.,., = ^ tr{VlV^,VlV^,) i- tr{VlV^,). (4.37) 

The choice is of direct phenomenological interest, since the operator above is the most compli- 
cated quantity appearing in the calculation of di-hadron production in proton-nucleus collisions^ 
[2-5]. Considering the same configuration as before, that is, taking ri3 = ri4 and r23 = ^24, we 
see that the evolution couples the operator in Eq. (4.37) to SxiX2SxiX2Sx2Xi- After a straight- 
forward calculation, similar to the one leading at Eq. (4.36), we arrive at 

{QxiX2XsX4Sx4X3)y = {SxiX2Sx3X4Sx4X3)y = {SxiX2)y {SxsX4)y'^ ~^ {SxiX2)y ■ 



(4.38) 



(Once again, we have assumed this equation to hold already in the initial condition at Iq, which 
is in particular true within the MV model, as one can check.) It is interesting that for this 
configuration, and similar to Eq. (4.36), the result depends only on ri2 and r34, but not on ri4 
and r23. Also, by keeping Nc finite, we notice that the last term in the r.h.s. (the single dipole 
S'-matrix) dominates deeply at saturation. This is in agreement with our earlier discussion in 
3.4, since the linear combination 

QX1X2X3X4SX4XS iV^ '5'a;ia;2 (4.39) 



In fact, the operator appearing in such a process is (5x20:3x4x1 <S'a;3a;4, but, due to the invariance under charge 
conjugation, its expectation value is equal to the one of the operator in Eq. (4.37). 
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is a special case of the operator in Eq. (3.37) for this particular configuration. In fact, the 
operator that appears in the di-hadron cross section is the combination 



- - 1 ' 

S(3X^X2XsXi ^^^^2 '^Qx\X2X^X4,Sx^X^ '^SxiX21 (4.40) 

which is also the quantity studied numerically in [11]. Using the previous results one finds 
that, for our special configuration, the expectation value of Sq is given by the relatively simple 
expression 

{3^x1x2x3x4) Y = {SxiX2)y {Sx-iXA)Y " ■ (4-41) 

It is amusing to note that for this particular configuration and for large iVc, the 4-point function 
relevant for di-hadron production factorizes as {Sqxix2X3X4)y = {Sxix2)y Sxzxi)^ ■ This is 
precisely the factorization formula used (for a generic configuration) in the phenomenological 
study in [6] — at that time, by lack of a better formula. Such a factorization however has no 
deep justification and is merely a property of the configuration at hand. As our next example 
will show, this 'factorization' can badly fail for other, equally simple, configurations. 

Specifically, let us consider the expectation value of the operator (4.37) for the 'line' config- 
uration studied in [11, 12] and shown in Fig. 3.d. Note that the two quarks of the quadrupole 
and the antiquark of the dipole are put in a same point {xi = x^), and similarly for the two 
antiquarks of the quadrupole and the quark in the dipole {x2 = x^). Thus, only one non-trivial 
distance r = ri2 = = ^34 = ri4 characterizes the configuration. Then one can easily check 
that the evolution of Qa;ia;2a;ia;2'S'a;2ici couples again to <S'a;ja;2'S'a;ia;2'^iC2a'i j leading to a 2x2 inho- 
mogeneous system of equations. Expressing 7y(r) in terms of the dipole (S'(r))y and using an 
obvious shorthand notation we have 



2 



- ^VI^T ^^^^^ + iV|3T is )y - {S)y, (4.42) 

'^''^^ - ' mY^'-§^iS^)Y-^iS)r. (4.43) 



d\n{S)Y N^--^' ' N^c-^ N^-^ 

The solution to this system is straightforward to obtain; so long as {QS)y is concerned, one 
finds 

(iVe + 2)(iV,-l) (iV, + l)(jV^-2) 

where we assumed that the above is already valid at Yq, as is the case in the MV model. Using 
this result together with Eq. (4.40), it is straightforward to evaluate {Sq)y for this particular 
configuration and compare with the numerical results in Ref. [11]. We shall find it rewarding to 
plot {Sq)y in two different ways; first as a function of rQg and then as a function of 1 — {S)y- To 
be in accordance with [11], the saturation momentum is defined by the condition {S)y = 
for rQs = a/2. 

We show this comparison in Fig. 4, where for some of the curves we have used the numerical 
data of [11]. On the left we show the correlator of interest as a function of rQs and for two 
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Figure 4. The expectation value of Sq, as defined in Eq. (4.40), for the 'hne' configuration. Left: as a 
function of the scahng variable rQs- Continuous magenta: JIMWLK at Y=0. Dashed cyan: MFA for 
Nc=3 at Y=0. Continuous red: JIMWLK at Y=5.18. Dashed blue: MFA for N^^S at Y^5.18. Right: 
as a function of 1— {S)y ■ Continuous red: JIMWLK at six different values of rapidity Y=0, 1.04, 2.07, 
3.11, 4.14 and 5.18. Continuous blue: complete result in the MFA for Nc = 3. Dashed magenta: large-Nc 
result in the MFA. Dotted gold: assuming factorization for the expectation value of QS and using the 
MFA for Nc = 3 for the expectation value of the quadrupole Q. Dotted dashed green: for illustrative 
purposes we also show {S)y. JIMWLK curves are constructed from the numerical solution given in [11]. 
MFA curves are analytical expressions in terms of (S), which is again provided by the numerical solution 
in [11] for the purposes of the left figure. 

values of the rapidity: y = 0, which is where one starts the evohition (with initial conditions 
of the MV type), and Y = 5.18, which is large enough for the effects of the evolution to be 
fully developed. The curves denoted as 'MFA' represent our present results, cf. Eq. (4.44) 
and (4.40), whereas the 'JIMWLK' curves follow from the numerical solution to the JIMWLK 
equation. For y = 0, the two types of curves overlap with each by construction, as they both 
reduce to the respective prediction of the MV model. What is remarkable though, is that a 
very good agreement between the numerical solution and the MFA persists for Y = 5.18. In the 
limiting regimes of weak and respectively strong scattering, the respectives curves are practically 
indistinguishable. In the transition region around rQs ~ 1, the agreement is not that perfect 
anymore, but the two curves are still very close to each other, confirming that the MFA is also 
an excellent global approximation. 

From Fig. 4 (left) we also notice that the shape of the curve changes as we evolve from 
y = to higher values of rapidity {Y = 5.18 in the figure). However, this change is mostly 
attributed to the evolution of {S)y as a function of rQs as the rapidity grows. Indeed, in the 
right panel of Fig. 4 we show the correlator of interest as a function of 1 — {S)y- One can see 
that the curves obtained from the numerical solution to JIMWLK for various values of Y form 
a very thin "band" whose borderline on the "lower" side is the MFA. This "band" is practically 
a line perfectly overlapping with the MFA when the scattering is either weak or strong, and 
becomes just a bit wider in the transition region. Thus, as we evolve in rapidity, the shape of 
the curve is barely changing. In fact, if one expands out the plot in order to better disentangle 
the various steps in the evolution, one can see that the high-y curve stabilizes very close to the 
y = curve after just a few units in rapidity. 
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Still in the right plot we also show two different approximations for {Sq)y, the one is the 
large-A^c result, while the other consists of factorizing {QS)y into {Q)y{S)y, as it would be 
justified at large Nc, but then using the finite- A'^c Gaussian approximation for {Q)y- This latter 
approximation takes into account some corrections, but not in a systematic way, and was 

used in [11] since the (MV-like) expression (4.44) was not available at the time. Comparing with 
the numerical findings in [11], we already saw that the complete result Eq. (4.44) at finite- A'^c 
is the one which shows the best agreement. Even though it is not very significant, we note 
that the large-A^c expression is the next one closer to the numerical data, perhaps because it 
is at least a systematic approximation. For illustrative purposes we also show (5')y, which is 
based neither on a large-A'^c approximation nor on a mean field one, but simply corresponds to 
a 'naive' counting of Wilson lines. It fails badly even in the BFKL regime and clearly it has no 
chance to describe properly the correlator of interest. 
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